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Abstract 

The  interactions  between  populations  can  be  positive,  neutral  or  negative.  Predation 
and  parasitism  are  both  relationships  where  one  species  benefits  from  t  he  interaction 
at  the  expense  of  the  other.  Predators  kill  their  prey  instantly  and  use  it  only  for 
food,  whereas  parasites  use  their  hosts  both  as  their  habitat  and  their  food.  I  am 
particularly  interested  in  microbial  parasites  (including  bacteria,  fungi,  viri,  and  some 
protozoans)  since  they  cause  many  infectious  diseases. 

This  thesis  considers  two  different  points  in  the  population- interaction  spectrum 
and  focuses  on  modeling  host-pathogen  and  predator- prey  interactions.  Tlie  first  part 
focuses  on  epidemiology,  i  e.,  the  dynamics  of  infectious  diseases,  and  the  estimation 
of  parameters  using  the  epidemiological  data  from  two  different  diseases,  phocine 
distemper  virus  that  affects  harbor  seals  in  Europe,  and  the  outbreak  of  HIV/ AIDS 
in  Cuba.  The  second  part,  analyzes  the  stability  of  the  predator-prey  populations 
that,  are  spatially  organized  into  discrete  units  or  patches.  Patches  are  connected  by 
dispersing  individuals  that  may,  or  may  not  differ  in  the  duration  of  their  trip.  This 
travel  time  is  incorporated  via  a  dispersal  delay  in  the  interpatch  migration  term,  and 
has  a  stabilizing  effect  on  predator-prey  dynamics. 

Thesis  supervisor:  Michael  G.  Neubert 

Title:  Associate  Scientist,  Woods  Hole  Oceanographic  Institution 
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Chapter  1 


Introduction 


1.1  Interactions  between  populations 

Populations  interact  in  different  ways.  Their  interactions  can  he  positive,  neutral 
or  negative.  Predation  and  parasitism  are  both  antagonistic  relationships  where  one 
species  benefits  from  the  interaction  at  the  expense  of  the  other.  Predators  kill  their 
prey  instantly  and  use  it.  only  for  food,  whereas  parasites  use  their  hosts  both  ns 
their  habitat  and  their  food.  Some  parasites,  including  bacteria,  fungi,  viri,  and 
some  protozoans,  are  microbial  and  replicate  within  the  host  (micToparasites).  Other 
pathogens,  such  as  ticks,  nematodes,  and  tapeworms,  have  no  within  host  replication 
and  are  known  as  macroparasites.  I  am  particularly  interested  in  microparasit.es  since 
they  are  pathogenic  agents  that  cause  many  infectious  diseases. 

Most  models  for  predator- prey  relationships  assume  prey  and  predators  have  rel¬ 
atively  equivalent  sizes  and  life  history  characteristics.  This  approach  is  applicable  to 
a  variety  of  organisms,  from  rabbits  and  foxes,  to  various  macroparasites.  But,  when 
the  host  is  very  large  with  relatively  slow  dynamics,  and  the  pathogen  is  small  and 
multiplies  rapidly  inside  the  host,  the  traditional  predator-prey  approach  is  not  very 
useful.  Here  the  dynamics  of  such  host- pathogen  interactions,  typical  for  infectious 
diseases,  is  studied  using  a  different  mathematical  approach.  Given  the  fast  dynamics 
of  the  pathogen,  and  the  slow  dynamics  of  the  host,  we  assume  that  the  host  pop¬ 
ulation  is  constant  and  we  ask  questions  about  the  spread  of  the  pal  hogen  between 
infected  and  not- infected  segments  of  the  host  population. 

In  this  thesis,  1  focus  on  two  extremes  along  this  consumer  resource  sprectrum 
and  contrast  host-pathogen  and  predator- prey  interactions. 
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Tlie  first,  part  (Chapters  2  -  4)  focuses  on  the  host-pathogen  end  of  the  spectrum 
with  Chapters  2  and  3  studying  the  Phocino  Distemper  Virus  that  affects  harbor 
seals  in  Europe,  and  Chapters  4  focusing  on  the  outbreak  of  HIV/AIDS  in  Cuba. 
The  second  part  {Chapters  5  and  6)  studies  theoretical  models  for  spatially-extended 
predator  prey  populations. 

In  this  introductory  chapter  I  discuss  the  basic  biology  of  Pliocine  Distemper 
Virus  outbreaks,  the  pathology  of  the  virus,  t  he  population  biology  of  harbor  seals, 
the  available  data  set,  and  the  motivation  for  the  spatial  predator  prey. 


1.2  Epidemiology 

An  outbreak  of  a  disease  that  spreads  rapidly  and  infects  a  substantial  portion  of 
the  population  in  a  region  over  a  short  period  of  time  is  known  as  an  epidemic. 
Major  epidemics  in  the  past  include  the  bubonic  plague  (“Black  Death")  that  spread 
from  Asia  throughout  Europe  in  14-th  century.  It  is  estimated  that  bubonic  plague 
has  killed  to  one  third  of  Europeans  between  1346  and  1350.  Another  example  is 
smallpox,  which  was  brought  to  North  America  by  invading  Spaniards  and  in  some 
cases  reduced  indigenous  population  to  one  tenth  of  it  s  preepidemic  size  t  he  Indian 
population  of  Mexico  is  thought  to  have  been  reduced  from  30  million  in  1510  to  only 
3  million  in  1530.  The  “Spanish  flu”  H1N1  influenza  epidemic  of  1918-1919  caused 
10-20  millions  of  deaths  worldwide. 

Major  epidemics  today  include  malaria,  tuberculosis  and  HIV/ A1  OS,  which  com¬ 
bined  kill  over  6  million  people  each  year.  At  the  end  of  2004,  there  were  40  million 
people  infected  with  HIV,  5  million  infected  in  2004  alone  (UNAIDS,  2004). 

When  ail  epidemic  affects  an  animal  population,  it  is  called  an  epizootic.  An 
example  of  a  recent  major  epizootics  is  Pliocine  Distemper  Virus  that  caused  the  death 
of  up  to  60%  of  harbor  seals  [Phoca  vitulina)  in  certain  locat  ions  in  the  Nort  h  Sea  in 
1988  and  2002  (R ducking.  2002,  2003;  Harkonen  cl  of.,  2006).  Other  species  within 
the  morbilliviridae  have  been  implicated  in  mass  mortalities  of  bottleuose  dolphins 
(Tursiaps  truncatus)  along  the  U.S.  Atlantic  coast  (1987-88)  and  Gulf  of  Mexico 
(1993-94),  striped  dolphins  (Stenella  coeruleoalba)  in  the  Mediterranean  Sea  (1990- 
92)  and  common  dolphins  ( Delphinus  dclphis)  in  the  Black  Sea  (1994)  (Osterhaus. 
1988:  Dietz  ct  ai,  1989:  Domingo  et  al..  1990:  Kennedy,  1998). 
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Figure  1-1:  Maps  of  1988  and  2002  PDV  outbreaks  in  northern  Europe.  Circles 
indicate  how  many  seals  died  in  a  particular  region.  After  Markon  et  al.  (2003). 


1.2.1  Phocine  distemper  virus 

Phocine  distemper  virus  (PDV)  was  first  described  in  1988.  when  it  killed  over  23, 000 
harbor  seals  (Phoca  vitulina)  in  northern  Europe  (Harkonen  et  al...  2006).  The  ‘seal 
plague  started  at  the  Danish  island  of  Anholt  (see  Figure  1-1)  anti  the  virus  quickly 
spread  to  populations  in  Sweden,  Netherlands,  England,  Scotland  and  Ireland  (Dietz 
et  al..  1989).  It  resulted  in  the  largest  recorded  epizootic  of  any  marine  mammal 
population  with  an  estimated  mortality  of  56-58%  in  large  regions  (Dietz  et.  al.  1989; 
Heide- Jorgensen  &  Harkonen,  1992;  Harkonen  et  al. ,  2002;  Harding  et  al..  2002). 

PDV  caused  another  outbreak  in  2002  (see  Figure  1-1),  killing  33.000  harbor 
seals  in  Baltic,  Wadden  and  North  Seas  between  May  and  October  (Reinekiug.  2002; 
Harkonen  et  al.,  2006).  Both  epizootics  originated  at  the  same  location,  the  island  of 
Anholt,  but  it  appeared  23  days  later  in  2002.  The  population  size  on  the  European 
continent  in  2002  was  about  twice  that  of  1988,  where  in  the  UK  it  was  at  comparable 
levels  both  years.  On  the  percent  basis,  tin1  two  outbreaks  had  comparable  mortality. 


1.2.2  Population  biology  of  harbor  seals 

The  dynamics  of  an  infectious  disease  is  determined  by  the  size  and  structure  of  the 
host  population,  its  life  history,  and  the  behavior  of  individuals.  All  of  these  factors 
affect  the  transmission  of  the  virus,  and  influence  the  dynamics  of  the  disease. 
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Figure  1-2:  Harbor  seal.  Pham  vitulina. 
http://www.wildlife.shetlaiul.co.uk/. 


Table  1.1:  Harbor  seal  population  sizes  before  and  after  the  1988  and  2002  epizootics 
in  Europe. 


Location 

Pre-1988 

Post-1988 

P  re-2002 

Post-2002 

Kattegat-Skagerrak 

12.700 

5.000 

23.000 

11,750 

Wadden  Sea 

10.840 

7,000 

35.000 

18.980 

European  continent 

30,800 

18.770 

02,070 

33.490 

UK 

53.000 

48,000 

53.000 

50,000 

Grand  total 

82.000 

59.500 

110.000 

82.800 

Population  size 


There  are  five  recognized  subspecies  of  harbor  seal  distributed  widely  over  the  north¬ 
ern  hemisphere.  Current  estimates  of  population  sizes  arc  imperfect  and  often  out¬ 
dated,  so  it  is  hard  to  say  how  many  harbor  seals  there  are  in  the  region.  In  the 
eastern  Atlantic  region  there  were  about  80.000  before  the  1988  epizootic,  of  which 
30.800  were  estimated  to  be  on  the  European  continent  (Harkoneu  v.t  al..  2000).  Ar¬ 
eas  of  greatest  abundances  were  Great  Britain  (53.000),  the  Wadden  Sea  (10,840) 
and  Kattegat-Skagerrak  (12.700)  (Harkoneu  ct  al..  2000).  After  the  PDV  outbreak, 
European  harbor  seal  population  dropped  to  an  estimated  59.500  (see  Table  1.1  for 
more  details).  During  the  next  14  years  the  population  recovered,  and  reached  a  size 
double  that  of  pre-1988  epizootic  on  the  continent  (02.000  in  2002),  and  comparable 
levels  to  the  pre-1988  epizootic  in  the  UK  (see  Table  1.1).  The  estimated  total  pop¬ 
ulation  size  in  2002  was  110,000  seals  (which  was  reduced  to  about  83.000  after  the 
second  PDV  outbreak). 
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Demography 

Common  seals  have  a  yearly  pupping  season  that  runs  from  late  May  to  early  August, 
and  peaks  in  late  June/early  July  in  Europe.  Females  mature  at  3  to  4  years;  males 
mature  a  year  later.  Mothers  suckle  their  pups  during  a  4-week  nursing  period,  after 
which  the  pups  undergo  a  post- weaning  fast  lasting  3  —  (i  weeks,  when  they  start  to 
catch  their  own  food.  After  the  pups  are  weaned,  seals  mate  in  water.  Fertilization 
is  followed  by  embryonic  diapause  (prolonged  period  of  delayed  implantation)  that 
husts  about  2.5  months.  The  total  gestation  period  is  around  10.5  months.  Pregnancy 
rates  exceed  809?  (Burns.  2002)  and  80  —  979?  of  the  mature  females  bear  a  pup  each 
year  (Harkbnen  k  Heide- Jorgensen,  1990).  The  single  pup  per  female  per  year  poses 
a  major  constraint  on  population  growth.  Heide-.Jorgensen  e.t  al.  (1992a)  constructed 
a  matrix  population  model  and  calculated  that  in  the  stable  age  distribution  the 
asymptotic  growth  rate  A  was  1.112.  Annual  survivorship  of  adult  harbor  seals  in 
the  absence  of  PDV  is  around  909?  for  males  (data  for  females  suggest  959?  survival) 
(Harkbnen  &  Heide- Jorgensen.  1990). 

Dispersal 

For  determining  epidemic  behavior,  besides  knowing  the  total  population  size,  it  is 
also  important  to  know  the  population  structure,  mixing  of  individuals,  as  well  as 
the  mixing  of  subpopulations.  Telemetry  studies  (Thompson  k  Harwood.  1990)  and 
long-term  study  of  freeze-branded  animals  (Harkbnen  k  Harding.  2001)  suggest  a 
high  degree  of  site  fidelity  among  adult  harbor  seals.  None  of  1G3  branded  seals  were 
observed  to  haul  out  beyond  a  32-km  radius  from  the  site  where  they  were  branded 
as  pups,  (Harkbnen  k  Harding.  2001).  On  the  smaller  scale,  there  is  a  strong  spatial 
segregation  by  age  and  sex,  as  well  as  different  migration  tendencies  between  sexes  and 
ages.  A  genetic  study  of  micro-satellite  polymorphism  (Goodman,  1998)  suggests  six 
distinct  population  units:  Ireland-Scotland,  English  east,  coast,  Wadden  Sea,  Western 
Scandinavia,  East  Baltic  and  Iceland. 

Haul-out  behavior 

Seals  give  birth,  rear  their  offspring  and  molt  on  land.  That  introduces  seasonality 
to  the  seal  haul-out  behavior,  which  peaks  during  the  pupping  and  molting  seasons 
in  late  spring  and  summer.  The  molt  occurs  during  midsummer  to  early  fall,  after 
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Cetaceans 


Figure  1-3:  Phylogenetic  tree 
showing  the  relationships  between 
the  different  morbilliviri  based  on 
part  ial  sequence  of  one  of  the  pro¬ 
tein  coding  genes  (the  P  gene). 
The  branch  lengths  are  propor¬ 
tional  to  the  mutational  differ¬ 
ences  between  the  viruses  and 
the  hypothetical  common  ancestor 
that  existed  at  the  nodes  in  the 
tree.  Source:  Barrett  (1999). 


cessation  of  the  breeding  season.  At  that  time  up  to  57 %  of  a  colony  can  be  found  on 
land  (Heide- Jorgensen  U  Harkoncn,  1992).  There  arc  differences  in  haul-out  timing 
among  age  and  sex  cohorts.  Yearlings  usually  molt  earliest,  followed  by  subadults, 
then  adult  females  and  adult  males  molt  last  (Burns.  2092).  Therefore,  the  number  of 
animals  on  land  depends  on  age,  sex.  and  the  time  of  year  (Thompson.  1939;  Harkoncn 
et  at.  1999.  2002).  Seals  also  haul  out  throughout  the  year  but  less  frequently  and  in 
smaller  numbers  than  in  the  summer. 

1.2.3  Phocine  distemper  pathology 

The  agent  responsible  for  mass  die-offs  of  seals,  phocine  distemper  virus  (Cosby  < i  at, 
1988;  Mahy  et  at.  1988:  Osterhaus  ct  at,  1989).  belongs  to  the  M.orbilliviridae  genus 
(see  Figure  1-3);  a  group  of  RNA  viruses  that  cause  infectious  diseases  in  mammals 
(Barrett  et  at,  1993;  Forsyth  et  at.  1998;  Barrett.  1999)  and  measles  in  humans 
(Barrett,  1987).  PDV  is  most  closely  related  to  the  canine  distemper  virus  (CDV) 
(Osterhaus  et  at.  1988;  Kennedy  et  at.  1988;  Rima  et  at.  1992)  that  can  cause  similar 
infections  in  other  seal  species  (Grachev  et  at.  1989;  Kennedy  et  at.  2000).  PDV  is 
thought  to  be  endemic  in  arctic  harp  seals  ( Phoca  groe.nlandica)  (Markussen  V  Have. 
1992).  In  1987  and  1988  harp  seals  were  observed  to  migrate  as  far  south  as  Danish 
waters,  and  are  thought  to  have  spread  the  virus  to  the  previously  unexposed  harbor 
seal  population  (Goodhart,  1988). 

The  disease  transmits  between  animals  in  close  contact  in  the  same  way  a  cold 
spreads  in  humans  -  by  inhalation  when  an  infected  individual  coughs  or  sneezes 
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(Kennedy.  1990.  1998).  Once  in  the  host,  the  virus  spreads  through  macrophages, 
lymphocytes  and  thrombocytes  and  infects  various  tissues.  Since  the  virus  is  spread 
in  air.  an  infected  animal  can  only  spread  the  disease  to  its  neighbors  during  haul 
outs. 

Symptoms  of  the  PDV  infection  in  seals  include  fever,  respiratory  problems  such  as 
coughing,  nasal  discharge,  as  well  as  discharge  from  the  eyes,  conjunctivitis.  Infected 
pregnant  females  abort  their  pups,  and  elevated  numbers  of  aborted  pups  in  a  certain 
area  can  be  an  indicator  of  the  presence  of  PDV.  A  pup  that  is  orphaned  or  abandoned 
before  weaning  will  die.  If  the  virus  enters  the  central  nervous  system,  infected  seals 
become  disoriented  and  will  be  disinclined  to  move  which  can  cause  them  to  spend 
more  time  on  land,  and  loss  time  in  water  searching  for  food.  Postmortem  findings 
inelm le  subcutaneous  emphysema  (air  bubbles  under  the  skin)  of  the  head  and  neck 
(Bergman  et  al..  1990;  Munro  et  al,  1992;  Baker.  1992)  due  to  which  dead  animals 
float  for  a  long  time  before  being  washed  up  on  land.  The  disease  is  confirmed 
by  blood  testing  diseased  animals  or  tissue  sampling  dead  animals  (Barrett.  1999). 
Morbilliviri  are  known  to  suppress  their  host's  immune  system,  thus  increasing  the  risk 
of  secondary  infection  by  a  wide  range  of  agents.  Autopsies  of  dead  seals  have  shown 
that  the  main  proximate  cause  of  death  is  a  secondary  infection,  bacterial  pneumonia 
caused  by  Dordetella  bjvnchiseptica  (Baker  k  Ross,  1992;  Kennedy,  1998), 

In  1988  the  disease  quickly  spread  from  central  Kattegat  (see  Figure  2-7).  where 
it  appeared  in  April,  to  Danish  and  Dutch  Wadden  Sea  (May),  and  then  to  Skagerrak 
and  German  Wadden  Sea  (June).  By  mid  July  seal  herds  in  the  Oslo  Fjord  and  the 
Baltic  were  affected.  British  liaul-oiits  were  the  last  to  be  bit  by  disease  in  August 
and  September  (Dietz  et  ai.  1989).  The  estimated  rate  of  spread  was  3.970  km/year 
(McCallum  et.  al.  2003).  In  each  location  the  epidemic  lasted  70-100  days;  longer 
when  haul-outs  were  less  discrete  as  in  the  Wadden  Sea  (Dietz  et  al.  1989) 


1.2.4  Data  on  haul-out  behavior 

Harbor  seals  have  been  st  udied  at  their  haul-out  sit  re  for  decades.  During  1978- 
1998  aerial  surveys  were  conducted  for  Swedish  and  Danish  haul-out  locations,  which 
were  photographed  in  the  peak  haul-out  season  and  seals  were  later  counted  from 
photographs  (Heide- Jorgensen  et  al.  1992b;  Harkonen  et  al..  1999,  2002).  The  sex- 
related  and  age-specific  seasonal  behavior  of  seals  has  been  inferred  by  studying  163 
freeze-branded  animals  during  1985-1997  (Harkonen  et  al. ,  1999)  and  8  VHF-taggcd 
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seals  during  1984-1986  (Thompson.  1989). 


1.2.5  Epizootic  data 

Major  haul-out  locations  in  the  Kat t egat-Skagerrak  were  regularly  surveyed  for  seal 
carcasses  by  biologists  and  other  trained  personnel  during  both  epizootic  periods 
(Harkdnen  k  Heide- Jorgensen.  1990).  In  the  UK.  most  of  tin1  dead  seals  were  re¬ 
ported  by  the  general  public  via  a  “hotline"  number  (http://www.defra.gov.uk/). 
The  number  of  reported  stranded  seals  are  treated  as  the  ‘number  of  dead  seals  per 
day'  which  form  the  epizootic  curves  (cumulative  numbers  of  dead  seals).  By  compar¬ 
ing  the  number  of  seals  found  dead  to  numbers  of  seals  in  population  surveys  before 
and  after  the  epizootic  we  can  estimate  the  mortality  in  each  location  (Dietz  ct  ai. 
1989;  Rcineking,  2002).  For  most  locations  the  day  that  first  dead  seal  appeared  is 
also  known. 

This  data-set  is  vulnerable  to  several  sources  of  observational  error.  The  number 
of  stranded  carcasses  depends  on  wind  directions  and  reporting  effort.  Carcasses  may 
float  for  a  while  before  finally  getting  washed  ashore,  so  the  seals  mav  be  lost  or 
washed  up  on  shores  far  away  from  their  actual  territory  and  dead  seals  might  be 
reported  several  times  (Thompson  k  Miller,  1992). 

The  epidemic  curves  are  the  only  epidemic  data  available  for  PDV  outbreaks,  and 
they  form  a  link  between  the  data  and  the  models. 

1.2.6  PDV  modeling 

Chapter  2  presents  a  model  for  PDV  outbreaks,  that  includes  the  information  oil  the 
life  history  of  seals,  and  the  transmission  of  the  virus.  Using  the  model,  I  develop 
a  way  to  estimate  epidemiological  parameters  based  on  the  available  epizootic  data. 
Seasonal  haul-out  behavior  influences  the  mixing  between  seals  and  the  transmission 
of  the  virus.  The  process  of  transmission  in  Chapter  3  incorporates  the  haul-out  be¬ 
havior,  and  I  use  this  model  to  investigate  differences  in  mortality  between  locations 


1.3  Prey-predator  interactions 

The  most  basic  models  (Lotka,  1926;  Vol  terra,  1931;  Nicholson  k  Bailey.  1935)  and 
experiments  (Cause.  1934)  predict  instability  of  predator  prey  systems.  How  do  then 
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predator  prey  systems  persist  stably  in  nature?  The  answer  most  often  given  is  that 
the  models  and  experiments  omit  processes  that  affect  stability  in  natural  systems 
(for  examples  see  May,  1973;  Hassell.  1978;  Crawley,  1992:  Mueller  k.  Joshi,  2000). 

Natural  systems  are  spatially  structured.  Populations  are  often  organized  into 
discrete  spatial  units  or  patches  that  arc  connected  by  dispersal  (metapopulation 
structure).  In  the  traditional  approach,  dispersal  is  assumed  to  occur  instantaneously, 
leaving  the  dynamics  of  the  model  often  unchanged.  In  reality,  individuals  spend  a 
finite  amount  of  time  in  transit  from  one  patch  to  another.  This  travel  time  can  be 
incorporated  in  predator  prey  models  via  a  delay  in  the  inter- patch  dispersal  term. 

What  are  the  effects  of  dispersal  delays  on  the  dynamics  of  the  predator  prey 
models?  To  find  out.  I  developed  predator  prey  models  that  include  dispersal  delays. 
These  models  have  tin1  form  of  a  system  of  delayed  differential  equations,  1  study  the 
dynamics  of  these  systems  analytically  and  numerically. 

To  determine  whether  the  dispersal  delays  have  a  stabilizing  effect  on  the  predator 
prey  equilibrium  point.  I  incorporated  dispersal  delays  into  the  Lotka  Volterra  model. 
The  equilibrium  point  of  the  non-spat ial  predator  prey  Lotka  Volterra  model  is  a 
center,  i.  e.,  a  "neutrally  stable’’  equilibrium  surrounded  by  a  family  of  periodic  solu¬ 
tions  whose  amplitudes  depend  on  the  initial  conditions.  The  slightest  change  to  the 
model’s  structure  typically  results  in  qualitatively  different  behavior.  For  example, 
if  the  growth  rate  of  prey  decreases  linearly  with  prey  density  the  equilibrium  point 
is  stable;  on  the  other  hand,  introducing  a  saturating  (Type  II)  functional  response 
turns  the  equilibrium  into  an  unstable  spiral  point  (Gotclli,  1995).  In  Chapter  5, 1  use 
this  structural  instability  of  the  Lotka  Volterra  model  to  show  that  dispersal  delays 
stabilize  the  equilibrium  point  of  the  spatially  structured  Lotka  Volterra  model. 

However,  the  the  Lot ka  Volterra  model  is  considered  oversimplified  for  two  rea¬ 
sons.  First,  in  the  absence  of  the  predators,  the  prey  grow  exponentially  without 
bound.  Second,  the  per  capita  rate  of  consumption  of  prey  by  predators  grows  in 
proportion  to  the  prey  population  size,  implying  that  individual  predators  can  process 
prey  items  infinitely  fast  .  These  faults  are  eliminated  in  the  Rosenzweig  MacArthur 
model,  which  includes  a  carrying  capacity  for  the  prey  and  a  finite  prey  handling  time 
for  the  predators  that  results  in  a  saturating  functional  response. 

The  Rosenzweig  MacArthur  model  has  a  more  complicated  dynamics  than  the 
Lotka  Volterra  model.  For  small  values  of  carrying  capacity,  the  coexistence  equilil>- 
riurn  point  is  locally  asymptotically  stable,  As  the  carrying  capacity  increases  beyond 
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some  threshold  value,  a  Hopf  bifurcation  occurs,  the  equilibrium  point  becomes  un¬ 
stable.  and  trajectories  are  drawn  onto  a  single  stable  limit  cycle.  This  destabilization 
by  increasing  prey  carrying  capacity  is  known  as  the  ’paradox  of  enrichment  ‘  (Rosen- 
zweig.  1971;  May,  1972:  Gilpin.  1972}. 

Dispersal  delays  are  strong  enough  to  overcome  the  destabilizing  effect  of  the 
Type  II  response  and  can  stabilize  the  coexistence  equilibrium  of  tin1  Rosenzweig 
MacArthur  model.  For  many  parameter  values,  stability  persists  even  in  the  limit  of 
infinite  carrying  capacity.  Dispersal  delays  also  help  resolve  the  paradox  of  enrichment 
by  reducing  the  amplitude  of  oscillations  when  the  equilibrium  is  unstable. 
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Chapter  2 


Estimation  of  the  basic 
reproductive  number,  TZq 


2.1  Introduction 

Phocine  distemper  virus  (PDV)  caused  mass  die-offs  of  European  harbor  seals  ( Phoca 
vitulina )  in  1988  and  2002  (Osterhaus  et  ai.  1988;  Maliy  et  ai,  1988;  Cosby  et  ai. 
1988;  Rima  et  at,  1992).  Botli  outbreaks  started  at  the  Danish  island  of  Anholt  in  the 
spring,  and  in  the  following  months  spread  throughout  the  entire  European  harbor 
seal  population,  killing  more  than  23,000  seals  in  1988  and  30.000  in  2002  (Harkonen 
ct  ai.  2000).  These  are  the  largest  epizootics  ever  reported  for  any  marine  mammal 
population  with  estimated  mortality  of  5G  -  58$?-  in  large  regions  (Dietz  ct  ai.  1989; 
Heide- Jorgensen  k  Harkonen.  1992;  Harkonen  et  ai,  2002;  Harding  ct  ai,  2002). 

It  is  likely  that  PDV  will  revisit  the  harbor  seals  of  Europe.  How  many  seals  will 
eventually  become  infected?  How  fast  will  the  epidemic  spread?  How  long  will  it  last? 
What  measures  should  be  taken  to  control  or  prevent  the  outbreak?  These  quantities 
are  related  to  the  epidemiological  parameter  known  as  the  basic  reproductive,  number. 
fto- 

Hit  is  defined  as  the  expected  number  of  new  infections  caused  by  a  single  infected 
individual  in  an  entirely  susceptible  population  (Dietz.  1975).  As  a  result.  Hu  provides 
a  threshold  for  whether  or  not  an  epidemic  will  occur.  If  Hq  >  1,  the  number 
of  infections  increases,  leading  to  an  epidemic;  if  Hq  <  1.  the  infection  dies  out 
(no  epidemic).  also  determines  the  duration  of  a  closed  epidemic  (Anderson. 
1990)  as  well  as  its  final  size  (Kermack  k  McKendrick,  1927;  Heesterbeek  k  Roberts. 
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1995;  Anderson,  199G).  The  quantity  77o  also  has  applications  in  developing  control 
strategies.  In  order  to  prevent  or  contain  an  epidemic,  the  proportion  of  individuals 
that  needs  to  be  removed  from  the  susceptible  pool,  either  by  vaccination  or  by 
culling,  is  1  -  l/770  (Anderson.  199G). 

The  concept  of  a  critical  threshold  arose  from  the  analysis  of  mathematical  models 
for  vector  borne  diseases  (Ross,  1911).  and  directly  transmitted  diseases  (Kern tack 
&  McKendrick,  1927)  at  the  beginning  of  the  last  century.  But.  it  was  not  until  the 
1980s  that  the  potential  of  use  of  77()  in  epidemiology  and  in  the  control  of  infec¬ 
tions  diseases  was  fully  recognized  (Dietz.  1975;  Diekmami  at  al..  1990;  Anderson. 
199G:  Dietz.  1993;  Heesterbeek,  2002).  This  is  surprising,  as  the  analogous  concept 
known  as  the  "net  reproductive  rate"  (denoted  by  77u  by  Dublin  k  Lotka  (1925)) 
was  fully  developed  for  the  study  of  demography  and  ecology  about  fifty  years  before 
its  widespread  application  in  epidemiology  (Sharpe  k  Lotka,  1911;  Dublin  k  Lotka, 
1925).  In  population  biology.  7Z{}  is  defined  as  the  expected  number  of  offspring  that 
an  individual  will  produce  during  its  lifetime,  or  the  population  growth  rate  from  one 
generation  to  the  next  (Caswell,  2001). 

In  tins  chapter.  I  focus  on  7 7q  in  two  ways;  (i)  given  a  model,  how  does  one 
calculate  77n .  and  (ii)  given  data,  can  77(]  be  estimated  directly,  or  how  does  one 
estimate  individual  parameters  in  the  model  to  determine  77().  In  order  to  estimate 
77|)  for  a  particular  disease,  one  needs  some  form  of  data  on  the  number  of  cases  that 
suffer  infection  from  this  disease.  Epidemic  data  from  naturally  occurring  wildlife 
diseases  is  often  lacking  in  detail  and  estimating  epidemic  parameters  is  challenging. 
In  the  ease  of  the  1988  and  2002  phocine  distemper  virus  outbreaks,  only  the  number 
of  stranded  seal  carcasses  were  observed.  Extensive  efforts  were  made  during  both 
the  1988  and  2002  epizootics  to  count  the  seals  that  died.  Time  series  of  stranded 
carcasses  collected  by  teams  of  biologists  and  other  trained  personnel  in  each  region 
were  used  to  construct  cumulative  curves  (also  called  epidemic  or  epizootic  curves). 
These  curves  form  a  link  between  the  data  and  the  models.  These  epizootic  curves, 
and  a  review  of  both  PDV  outbreaks  can  be  found  in  Harkonen  at  al.  (2006). 

In  the  next  section  I  present  a  short  review  of  different  calculations  of  77,, ,  and 
of  common  ways  to  estimate  77n  from  data.  Based  on  a  model  for  PDV  dynamics, 
I  develop  a  new  likelihood-based  method  for  estimating  77()  from  epidemic  curves.  1 
will  use  simulation  results  to  evaluate  accuracy,  precision,  and  the  bias  of  I  lie  method 
Using  this  method.  I  estimate  77()  values  for  different  regions,  and  show  that  regional 
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differences  in  71q  are  significant.  Further,  I  investigate  the  relationship  of  TZq  with 
variables  that  most  commonly  influence  it,  such  as  population  size,  spatial  structure, 
timing  of  infection,  and  the  level  of  immunity. 


2.2  Calculation  of  7Zu 


Consider  a  simple  deterministic  model  where  a  population  of  size  N  is  divided  into 
three  epidemic  compartments:  susceptible  individuals  S.  infective  individuals  /  (i. 
individuals  that  are  infectious),  and  a  removed  class  it  that  consists  of  individuals 
that  were  infected  but  are  no  longer  infectious  or  susceptible  to  reinfection  Such 
models  arc  often  called  SIR  models.  Contacts  between  individuals  are  assumed  to 
be  made  at  random.  The  disease  spreads  when  an  infectious  individual  contacts  and 
infects  a  susceptible.  The  force  of  infection  A.  defined  as  the  probability  per  unit  time 
for  a  susceptible  to  become  infected  (Diekmann  k  Heesterbeek,  2000),  is  a  product 
of  three  parameters:  (i)  tin1  contact  rate  c(N).  (ii)  the  probability  that  the  contact 
is  with  an  infective,  usually  assumed  to  be  f/N ,  and  (iii)  the  probability  p  that  the 
contact  between  susceptible  and  infective  individuals  in  fact  leads  to  transmission  of 
the  pathogen,  i.  e.,  the  probability  that  the  contact  is  'successful', 


A  .  .  c(JV)/p 


(2.1) 


The  per  contact  probability  of  successful  transmission  p  is  usually  assumed  to  be 
constant. 

If  the  rate  of  contacts  for  a  given  susceptible  individual  is  proportional  to  the 
population  size  N .  c(N)  —  C\N,  the  force  of  infection  is  proportional  to  the  number 
of  infect ives  l 

X~Cipf  =  j3I.  (2.2) 

The  proportionality  constant.  ft,  is  called  the  transmission  rate;  it  consists  of  both 
tin1  contact  rate  and  the  probability  that  the  contact  is  successful.  The  rate  at  which 
new  infections  occur  is  assumed  proportional  to  the  number  of  susceptibles  S.  and  is 
given  by  the  product  A S  =  j3SI. 

If  individuals  recover  from  infection  and  become  immune  at  a  rate  7.  we  obtain 
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the  standard  SIR  model 


^  =  -&SI,  (2.3a) 

at 

%  =  PS1  —  7/,  (2.3  b) 

%  =  7/  (2.3c) 

a/ 

first  studied  by  Kermack  k  MeKendrick  (1927).  (A  detailed  analysis  of  this  model 
can  be  found  in,  e.  <)..  Heesterbeek  k  Roberts  (1995);  Diekimum  k  Heesterbeek  (2000) 
or  Brauer  it  Castillo-Chavez  (2001).) 

Under  model  (2.3)  a  typical  infective  individual  meets  and  infects  .IS  susceptible 
individuals  per  unit  time,  and  continues  to  do  so  during  its  expected  infective  period 
1/7.  so  tin1  total  number  of  secondary  infections  that  individual  produces  is  ..IS’/ 7. 
If  at  the  beginning  of  the  epidemic  there  are  N  susceptible  individuals,  the  basic 
reproductive  number  for  model  (2.3)  is 

3N 

Tlo  =  — ■  (2.4) 

7 

This  type  of  transmission,  where  the  rate  of  contacts  for  a  given  susceptible  is  pro¬ 
portional  to  the  population  size  A,r.  is  known  as  mass  action  or  density- dcpendtui 
transmission  (Begun  id  al...  2002;  Brauer.  2006). 

When  the  number  of  contacts  per  infective  per  unit  time  is  constant.  c(,Y)  =  c2. 
the  process  of  transmission  is  known  as  pseudo-moss- action  (<■:.  3. in  Swinton  c.t  al.. 
1998).  standard  incidence  (Brauer,  2006).  or  fjcque.ncy-dcpcndt  nt  transmission  (Thrall 
k  Antonovics,  1997;  Begun  et  at.  2002).  This  type  of  transmission  is  most  commonly 
used  in  modeling  sexually  transmitted  diseases.  The  force  of  infection  for  frequency- 
dependent  transmission  is 

a  T  §■  (2-5) 

where  J  is.  again,  the  transmission  rate,  but  has  different  units  than  in  equation  (2.2). 
In  this  case,  the  basic  reproductive  number 

7Zo  =  -  (2.6) 

7 


is  independent  of  population  size. 


2.2.1  Survival  function 


A  more  general  formulation  of  1Z 0  follows  directly  from  its  definition.  Let  1(a)  be  the 
probability  that  a  newly  infected  individual  remains  infectious  for  at  least  time  a, 
and  let  m(a)  be  the  rate  of  infectious  ness  by  an  individual  that  has  been  infectious 
for  a  units  of  time.  The  number  of  secondary  infections  is  then  given  by 


(2.7) 


An  identical  formulation  is  found  in  population  biology,  where  1(a)  is  survivorship, 
t hat  is,  the  probability  of  surviving  from  birth  to  age  a.  and  ;ri(n)  is  the  rate  of 
reproduction  at  age  a  (e..g.  Heesterbeek  k  Roberts,  1995;  Keeling  k  Grenfell.  2000; 
Caswell,  2001). 

For  model  (2.3),  the  duration  of  infection  is  exponentially  distributed  with  the 
mean  1/7.  so  1(a)  =  exp(— 7a).  As  the  transmission  rate,  /?,  does  not  depend  on 
how  long  individuals  have  been  infectious,  the  rati'  of  infection  is  the  same  for  all 
infectious  individuals.  Overall  rate  of  infection  in  (2.3)  is  ,3SI ,  or  m(a)  =  liS  per 
infective  individual,  i.  e.,  /3S0  at  the  beginning  of  the  infection.  Substituting  for  1(a) 
in  m(a)  in  equation  (2.7)  gives 


(2.8) 


This  mathematically  natural  definition  of  7^  is  not  always  useful  for  computations, 
especially  when  dealing  with  more  complex  models 

2.2.2  Next  generation  method 

In  many  cases  it  is  useful  to  distinguish  between  different  classes  of  infectives  For 
example,  in  the  model  for  the  HIV/AIDS  epidemic  in  Chapter  4.  there  are  three  infec¬ 
tive  compartments  undiagnosed  HIV  cases,  diagnosed  HIV  cases,  and  AIDS  cases. 
In  other  situations,  multiple  infective  classes  can  be  used  to  capture  the  underlying 
age  structure  or  spatial  structure.  For  this  type  of  model.  can  be  derived  using 
the  next  generation  method  (Diekmann  et  a!..  1990;  de  Jong  et  ai„  1994;  Diekniami 
k  Heesterbeek,  2000;  van  den  Driessche  k  Watmough,  2002),  where  7?u  is  given  by 
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the  spectral  radius,  p,  (dominant  eigenvalue)  of  the  next  generation  matrix.  FV  1 


Ko=fi[FV-1]. 


(2.9) 


There  is  an  analogous  expression  in  discrete  time  demographic  models  of  the  form 
n(l  +  1)  =  An(;()  u(l),  The  vector  n(t)  describes  the  state'  of  the  population  at  time  /, 
and  let  the  projection  matrix  An  consist  of  the  transition  matrix  T„  and  reproduction 
matrix  Fn,  so  that  An  =  Tn  +  Fn.  Then 


Ha  =  p  [Fn(I  -  Tn)_l]  . 


(2.10) 


To  find  the  next  generation  matrix  FV~K  first  assume  there  are  n  compartments 
of  which  m  are  infective.  Lot  x  =  ,xn  f)e  the  number  of  individuals  in  each 

compartment.  Let  Ti{x)  be  the  rate  at  which  newly  infected  individuals  enter  com¬ 
partment  /.  V/  (x)  be  tin1  rate  of  transfer  of  individuals  into  compartment  /  by  all  other 
means  (including  the  transfer  of  infectious  individuals  from  one  infective  compartment 
to  another),  and  V~(x)  be  the  rate  at  which  individuals  arc  leaving  compartment  i. 
Define  V,(.r)  as  =  V“(z)  —  V,'  (x).  Tin-  rate  of  change  of  compartment  i  is  then 

X)  =  Ti  —  V,(.r).  We  can  then  form  the  next  generation  matrix  F I*"1  by 


(2.11) 


where  Lj  =  1 . m  and  j:()  is  the  disease  free  equilibrium,  at  which  the  popula¬ 

tion  remains  in  the  absence  of  the  disease.  (A  detailed  description  of  assumptions, 
constraints  and  proofs  of  theorems  can  be  found  in  van  den  Driessehe  &  Watmough 
(2002)).  The  (j.  k)  entry  of  V~]  is  the  average  amount  of  time  an  infective  individual 
that  was  introduced  into  compartment  k  spends  in  compartment  j  during  its  lifetime. 
The  (Lj)  entry  of  F  is  the  rate  at  which  infected  individuals  in  compartment  j  pro¬ 
duce  new  infectious  in  compart  merit  i.  Therefore,  the  entry  ( i.k )  in  the  generation 
matrix  /'V-1  is  the  expected  number  of  new  infections  in  compartment  /  produced 
by  an  individual  originally  introduced  in  compartment  k 

To  illustrate  this  approach,  imagine  a  case  where  the  early  stage  and  late  stage  of 
infection  have  different  transmission  rates,  0i  and  Lh-  We  can  model  this  scenario  by 
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a  simple  SIR  model  (or  S1IR).  with  two  infectious  compartments  /(  and  /jp 


J  (T 

—  ~(0\  h  + 

=  {fiyh+Mzh)? — 

dh  }  , 

ir  =  -i'1' 

dR 


dl 


=  72 1-2  ■ 


(2.12a) 

(2,12b) 

(2.12c) 

(2.12d) 


All  newly  infected  individuals  enter  the  compartment  l\  so  T\  =  (Ti  /(  +  .i2lz)S.  and 
?i  —  0.  Other  movements  among  the  compartments  are  described  by  Vi  =  7j  /]  and 
V2  =  72/2-  7iA,  giving 


F  = 

1 

£ 

i  - _ 

and  V  = 

7i 

0 

0  0 

.  “7! 

72  _ 

FU"1  - 


and  -  50  (ff  +  $)■ 


fli  ^'o  ^0 

7i  72  7a 


(2.13) 

(2.14) 


2.3  Estimation  of  TZu  from  data 

Contact  rates  and  transmission  rates  are  often  difficult  to  determine  from  observa¬ 
tions.  As  a  result,  7?0  is  difficult  to  calculate  using  equations  (2.7)  or  (2.9).  In  this 
section,  I  review  some  alternative  approaches  to  estimating  F0  from  data.  These  ap¬ 
proaches  either  assume  an  epidemic  in  a  closed  population  where  the  infection  leads 
to  immunity  or  death,  or  an  endemic  equilibrium. 


2.3.1  Final  size  equation 

In  the  case  of  a  closed  epidemic,  there  is  110  influx  of  susceptible  hosts.  Unlike  the 
endemic  ease,  where  a  pathogen  becomes  established  in  a  host  population,  in  a  closed 
epidemic  the  number  of  infections,  after  an  initial  increase,  eventually  drops  to  zero. 
The  fraction  of  the  individuals  that  eventually  become  infected  during  the  epidemic, 
the  final  size  of  the  epidemic,  was  first  analytically  determined  by  Kerinaek  k  MeK- 
endrick  (1927).  For  model  (2.3)  we  can  calculate  the  final  size  by  formally  dividing 


35 


equation  (2.3b)  by  equation  (2.3a),  and  integrating  the  expression  for  dl/dS  to  find 
the  orbits  in  the  (.S’,  /)  plane 


dl  =  ^-1  +  j  dS,  (2.15a) 

I  =  -S  +  ?  In  S  +  c,  (2.151)) 

where  c  is  the  arbitrary  constant  of  the  integration.  In  other  words,  orbits  are  defined 
by 

V\ S,  l)  =  S  +  I  -  -  hi  S  =  c.  (2 . 1  (i) 

0 


Since  none  of  these  orbits  reaches  the  I- axis.  S(t)  >  0  for  all  times.  The  part 
of  the  population  that  escapes  infection  is  Sl3C  =  liint_oc  5(f).  At  the  beginning  of 
the  epidemic  (t  =  (1).  all  of  the  individuals  are  susceptible  (S(0)  --  A'),  and  there 
an1  essentially  no  infected  individuals  {/(())  ~  0).  As  the  disease  disappears  from 
the  population  after  some  time,  there  are  no  infectious  individuals  at  the  end  of  the 
epidemic,  so  =  li*nt_oc  1(1)  =  0.  The  relation  V(Sq,  lQ)  —  V(SX.  Ix)  gives 


So  -  In  Sq  =  Sx  -  -  In  Sx, 

0  U 


In  —  =  -S0 


foc 

In  ^  =  1Z< 


S* 

So 


Sn 


=  exp 


7£d 


(2.17a) 

(2.17b) 

(2.17c) 

(2.17d) 


The  final  size  of  the  epidemic.  /.  is  simply  /  -  1  —  S^/So  —  X— 3(oo).  By  rearranging 
equation  (2.17e).  we  can  obtain  a  relationship  for  TZU  and  /. 


=  1»(1  -  /) 


(2.18) 


For  large  values  of  TZU. 


f  fa  l-e“%'  (2.19) 

is  a  useful  approximation  of  the  final  size  (Figure  2-2).  For  small  TZ{,.  the  final  size  is 
approximated  by 


j  ~  2(7^o  —  1 ) 


(2.20) 


3G 


5 


Figure  2-1:  Estimates  of  7£(1  using  the  final  size  equation  (2,18)  for  different  final  sizes 
of  the  epidemic. 

whit'll  follows  from  the  Taylor  expansion  of  (2.17e)  around  =  1  (or  s(oc)  =  1) 
(Diokmann  &  Heesterbeek.  2000.  pi 82). 

The  use  of  the  final  size  equation  (2.18)  for  estimating  7l0  requires  knowledge  of  the 
final  size  of  the  epidemic.  The  proportion  of  the  population  that  is  exposed  to  i  lie  virus 
during  the  outbreak  can  be  empirically  determined  by  extensive  serological  studies, 
For  the  PDV  outbreak  in  harbor  seals  in  Europe  such  studies  are  rare,  although 
some  morbillivirus  antibody  prevalence  studies  were  done  on  Scottish  populations 
(e.g.  Thompson  at  at. ,  1992,  2002;  Pomeroy  at  at .  2005),  suggesting  that  the  final 
size  of  the  1988  PDV  outbreak  in  Scotland  was  0.5-0. 7. 

Figure  2-1  shows  illustrates  the  relationship  given  in  equation  (2.18)  over  the  range 
of  all  possible  final  sizes.  Note  that  for  /  near  1.  small  errors  in  the  estimate  of  /  will 
produce  large  errors  in  the  estimate  or  TZU 

2.3.2  Proportion  of  susceptibles  at  the  endemic  equilibrium 

When  a  pathogen  invades  a  host  population,  the  number  of  susceptible  hosts  de¬ 
creases.  Eventually  the  system  may  arrive  at  an  equilibrium,  where  the  rate  at  which 
susceptible  individuals  are  infected  is  exactly  balanced  by  the  rate  at  which  newly 
susceptible  hosts  enter  the  population  (either  by  birth,  immigration,  or  loss  of  immu¬ 
nity).  This  is  known  as  an  endemic  equilibrium,  where  the  number  of  susceptible  and 
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Figure  2-2:  Comparison  of  the  true  final  size  and  its  approximation  given  with  equa¬ 
tion  (2.19). 

infective  individuals  is  given  by  (.s.  i)  —  (s*.  /*).  As  the  population  dot’s  not  consist  en¬ 
tirely  of  susceptible  individuals.  7?(,  is  not  directly  applicable  to  this  case.  Instead,  we 
can  use  the  effective  reproductive  number,  It,  which  is  the  number  of  secondary  cases 
produced  in  a  population  not  consisting  entirely  of  susceptible  individuals  (Grenfell 
&:  Dobson.  1995).  At  equilibrium,  each  infection  will  on  average  result  in  exactly  one 
secondary  infection,  so  the  effective  reproductive  number.  It.  will  be  equal  to  unity. 
In  a  homogeneously  mixed  population,  the  number  of  secondary  infect  ions  will  be 
proportional  to  the  probability  that  an  infectious  individual  contacts  a  susceptible. 
In  this  case,  the  effective  number  It  is  the  basic  reproductive  number  decreased 
by  a  fraction  of  host  population  that  is  still  susceptible,  ,s,  It  =  7Z{ts  (Dietz.  1975: 
Anderson,  1996).  At  equilibrium,  It  =  1  and  s  =  $*,  leading  to 


1 


(2.21) 


3* 


2.3.3  Average  age  at  infection 

Fur  a  homogeneously  mixed,  stationary  population  (where  births  exactly  balance 
deaths),  at  an  endemic  equilibrium,  Kq  can  be  estimated  as 


(2,22) 
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Table  2.1:  Summary  of  TZ^  calculations. 

Tii)  =  description  equation  number 


JSuh 

standard  SIR  model 

(2.4) 

Jo*  l{a)  in  (a) da 

survival  function 

(2.7) 

p{FV~l) 

next  generation  method 

(29) 

/ 

final  size  equation 

(2.18) 

1/s* 

endemic  equilibrium 

(2.21) 

L/A 

average  age  at  infection 

(2.22) 

Here  L  is  life  expectancy,  A  is  the  mean  age  at  infection,  and  the  infection  is  immuniz¬ 
ing  {Dietz,  1975:  Anderson,  1996).  If  the  net  population  growth  rate  is  positive,  the 
use  of  equation  (2.22)  can  lead  to  overestimation  of  7 7y.  as  the  relevant  demographic1 
time-scale,  i,  r...  the  reciprocal  of  the  birth  rate,  is  shorter  than  life  expectancy,  espe¬ 
cially  when  the  birth  rate  is  high.  For  endemic  infectious  in  a  growing  population,  a 
better  approximation  is 

(2.23) 

where  H  is  the  reciprocal  of  the  average  birth  rate  (Anderson.  1996). 

For  all  approximations  in  (2 . 18 )-{2. 23) .  one  must  also  assume  homogeneous  mixing 
in  the  population.  Therefore,  estimates  of  77o  obtained  from  (2.21)  or  (2.22)  would 
not  be  reliable  in  the  case  where  heterogeneity  in  host  demography  affects  the  contact, 
process  and  the  transmission  of  infection. 

Phocine  distemper  virus  is  not  endemic  in  harbor  seals  the  infection  seemed 
to  disappear  from  the  seal  population  following  both  the  1983  and  2(K)2  outbreaks 
so  we  cannot  use  equations  (2.21 )  (2.23)  to  estimate  7Zt}  PDV  outbreaks  are 
short  relative  to  the  life-span  of  this  long-lived  species,  so  we  can  assume  that  the 
population  did  not  grow  during  the  outbreak  and  consider  both  outbreaks  as  closed 
epidemics.  In  theory,  we  can  use  the  final  size  equation  (2.18)  to  roughly  estimate 
TZq.  In  practice,  the  data  on  the  fraction  of  the  population  that  escapes  the  infection 
comes  from  serological  surveys.  Such  surveys  are  rare  for  European  harbor  seals, 
so  there  is  not  enough  data  to  estimate  and  compare  7Zq  for  phocine  distemper  for 
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different  locations  in  Europe.  Approximations  of  7£0  summarized  in  Table  2.1  do  not. 
apply  for  PDV,  so  we  need  to  develop  a  new  method  to  estimate  Tin  from  the  data 
that  is  available  the  time-series  of  the  number  of  seals  that  have  died  from  PDV. 


2.4  Phocine  distemper  virus  dynamics 


The  methods  mentioned  in  Section  2.3  assume  that  the  epidemic  process  is  deter¬ 
ministic.  In  reality,  there  are  many  random  perturbations  that  can  influence  the 
transmission  of  the  disease  and  the  final  size,  and  those  random  effects  can  be  par¬ 
ticularly  important  for  small  compart mental  sizes,  or.  in  the  ease  of  seals,  for  small 
haul-out  units.  In  those  eases,  stochastic  models  are  more  appropriate.  Stochastic,  or 
probabilistic,  models  allow  for  the  use  of  more  sophisticated  estimation  procedures, 
such  as  maximum-likelihood-type  methods  that  require  data  on  the  time-series  of  the 
number  of  individuals  in  epidemic  compartments. 

The  time-scale  of  phocine  distemper  virus  outbreaks  is  much  shorter  than  the  de¬ 
mographic  time-scale  of  harbor  seals.  Therefore,  we  can  assume  that  seal  populations 
do  not  grow  during  an  outbreak,  and  we  approximate  a  single  PDV  outbreak  in  a 
single  haul-out  location  as  a  closed  epidemic.  We  take  a  compart  mental  approach  to 
modeling  PDV  dynamics,  and  we  divide  tin1  population  on  day  I  into  susceptible  seals 
($t),  infectious  seals  (/,}.  and  a  removed  class  [fit)-  which  consists  of  both  immune 
and  dead  seals. 

Let  the  number  of  seals  that  become  infected  on  day  I  (the  incidence)  be  a  random 
variable  A).  Let  xj  be  the  realization  of  the  random  variable  Xt.  i.c.  the  actual 
number  of  seals  infected  on  that  day.  After  a  seal  lias  been  infected  with  PDV.  we 
assume  it  experiences  a  latent  period  of  3  days  (Osterhaus  c.t  ai,  1938,  1989c;  Harder 
et  ai.  1990).  during  which  it  is  not  infectious.  Individuals  going  through  the  latent 
period  constitute  the  exposed  class,  which  is  not  directly  modeled  in  (2.24).  The 
model  accounts  for  the  exposed  class  by  introducing  a  delay  the  newly  infect ed 
individuals  enter  the  infectious  class  three  days  after  they  got  infected.  The  infectious 
period  lasts  12  days  (Osterhaus  et  ai.  1989a, c;  Grachev  c.t  ai.  1989;  Harder  c.t  ai. 
1990).  after  which  a  seal  cither  becomes  immune  or  dies.  Thus,  for  given  initial 
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conditions  we  can  compute  the  epidemic  trajectory  according  to 


Si+ 1  =  at-Xt, 

h+l  —  H  +  ^f-3  — 
R-t  +  1  —  rt  +  ]  5  > 


(2.24a) 

(2.24b) 

(2.24c) 


where  .r,  is  zero  for  /  negative  (Heide- Jorgensen  Harkdnen,  1992). 

An  individual  that  is  susceptible  at  time  I  remains  susceptible  at  time  I  +  1  only 
if  it  avoids  infectious  contact  with  all  it  infectives.  Let  p  be  the  probability  that  a 
given  infectious  seal  infects  a  given  susceptible  during  one  day,  The  probability  that 
a  susceptible  does  not  get  infected  upon  contact  with  a  given  infective  is  then  1  —  ]>. 
hence  (1  -  p)H  is  the  probability  of  not  getting  infected  by  any  of  the  i,  infectives 
at  time  /.  The  total  probability  that  a  susceptible  gets  infected  on  day  /  is  then 
1  -  (1  —  pY' .  and  this  event  is  independent  for  each  of  the  st  susceptible  individuals. 
Thus,  the  random  variable  X,  that  describes  new  infections  is  binomially  distributed. 


A',  ~  Bin[stl  1  -  (1  -  p)*‘j. 


(2.25) 


This  type  of  formulation  dates  back  to  the  series  of  lectures  by  Reed  and  Frost  in 
1928  (first  published  by  Abbey.  1952).  The  probability  of  having  an  epidemic  will  be 
a  product  of  a  chain  of  binomial  probabilities  of  the  form  (2.25).  Hence,  this  type  of 
a  model  is  referred  to  as  a  chain  binomial  model  (e.g.  Bailey,  1957;  Daley  Gaui. 
1999;  Andersson  Britton,  2000). 

We  further  assume  that  the  number  of  seals  that  die  each  day,  V).  also  is  binomially 
distributed,  with  constant  probability  of  dying  m. 


Yt  ^  Bin[xt_t5,m]. 


(2.20) 


For  stochastic  epidemics  in  large  communities  one  of  two  scenarios  can  occur. 
Either  a  small  number  of  individuals  get  infected,  or  there  is  a  major  outbreak. 
If  7lt)  <  1.  a  small  outbreak  may  occur;  major  outbreaks  arc  possible  if  and  only 
if  TZq  >  1  (as  stated  by  the  threshold  theorem  Bartlett,  1960).  In  this  case,  the 
asymptotic  distribution  of  final  sizes  of  the  epidemic  consists  of  two  parts,  first  one 
close  to  zero,  and  the  second  one  spread  around  the  deterministic  final  size  value  (see 
histogram  in  Figure  2-3). 
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Figure  2-3:  Simulated  epidemic  trajectories  and  histograms  of  final  sizes  of  1000 
simulations,  for  initial  susceptible  population  Sq  —  1000.  iQ  —  1.  m  —  O.G.  A)  and  B) 
TZ{)  =  0.8  C)  and  D)  1 Zq  -  1.5;  E)  and  Fj  TZq  =  3.  Dotted  lines  represent  20%  of  the 
expected  deterministic  final  size  approximated  by  (2.10). 
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2.4.1  Pseudo-maximum-likelihood  method  for  estimating  TZu 

Fur  model  (2,24)  the  basic  reproductive  number  is  given  by 


Ho  =  12pso,  (2.27) 

where  So  is  the  initial  susceptible  population  size.  Estimating  Ho  amounts  to  esti¬ 
mating  the  parameter  p. 

The  probability  of  infection  p  could  be  estimated  by  maximum  likelihood  methods, 
if  we  knew  the  number  of  seals  in  each  compartment  for  every  day  of  the  epizootic.  To 
see  this,  let  I  =  0  be  the  beginning  of  the  disease  outbreak,  and  let  /  =  T  indicate  the 
day  the  outbreak  ends.  Had  we  observed  the  number  of  seals  in  each  class  throughout 
the  outbreak,  the  likelihood  of  the  observed  trajectory  would  be 

r 

%)= (2-28} 
k-  0 

where  /  is  the  binomial  probability  density  function. 

f(x\s,p)=  -jp)*-*  (2.29) 

The  maximum  likelihood  estimate  of  p  would  then  be  the  value  p  which  maximizes 
Up). 

However,  as  in  most  wildlife  disease  outbreaks,  the  numbers  of  susceptibles  and 
newly  infected  etc.  were  not  observed  every  day.  Our  only  observation  is  the  number 
of  stranded  dead  seals  Each  day  a  certain  number  of  seals  dies  from  PDV  Out  of 
this  total  number  of  victims,  a  certain  proportion  strands  ashore,  and.  finally,  some 
proportion  of  tire  stranded  carcasses  gets  reported.  The  daily  number  of  reported 
stranded  carcasses  is  the  only  available  information  on  the  number  of  seals  that  died 
each  day  The  number  of  seals  that  gets  stranded  and  reported  will  depend  on  many 
factors  such  as  the  weather  conditions,  direction  of  the  wind,  stranding  location, 
reporting  effort,  etc.  Therefore,  there  is  inevitably  an  observational  error  associated 
with  the  daily  number  of  reported  carcasses,  and  the  observational  error  will  vary 
from  day  to  day. 

I  account  for  the  observation  error  in  a  simple  way.  There  are  two  sources  for  the 
information  on  the  total  number  of  seals  that  have  died  in  the  outbreak.  One  source 
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is  the  sum  of  the  daily  number  of  reported  strandings.  Themore  reliable  data  on  the 
total  number  of  seals  that  died  comes  from  census  data.  Since  we  know  the  total 
population  size  before  and  after  the  epizootic,  we  can  use  the  difference  in  census 
data  to  infer  the  total  death  toll.  I  model  the  observation  error  as  equal  to  the 
ratio  of  the  total  number  recovered  stranded  seals  to  the  difference  in  census  data.  1 
further  assume  that  the  observation  error  is  constant  throughout  the  outbreak  and 
scale  the  epidemic  curve  to  match  the  total  number  of  seals  that  died  according  to  the 
population  counts.  The  scaled  daily  counts  of  dead  seals  are  then  used  to  (Obstruct 
tin*  estimates  of  st,Xt  and  it. 

To  estimate  the  series  of  incidence,  we  equated  the  observed  mortality  with  its 
expectation  under  (2,20)  and  find 


m 

The  incidence  eases  can  only  have  integer  values.  Rounding  of  the  series  to  the 
nearest  integer  introduces  the  possibility  that  the  number  of  the  total  individuals 
infected  is  larger  than  the  initial  susceptible  population  size.  Therefore,  to  keep  .f 
series  in  integer  form.  I  round  the  right-hand  side  of  the  equation  (2.30)  to  the  nearest 
integers  towards  minus  infinity  using  the  Matlab  command  floorQ. 


With  j'i  in  hand  we  used  model  (2.24)  to  reconstruct  the  series  for  susceptible 
seals 

St+i  =  .s,  =  Af.  (2.31) 

and  the  series  of  infectious  seals  via 


it  —  it  1 1 


Vi  4i  2  —  ih  -  n 

— - — - — it  =  o 


(2.32) 


m 

We  then  treat  estimates  s ,x.  i  as  though  they  were  actual  observations  and  estimated 
p  by  minimizing  the  negative  log-likelihood  of  the  estimates 

r 


f(p)  =  -  ^2  1»  f  (xkrh,  1  —  (1  — 


fc-0 

over  p.  Tilt*  estimate  for  the  basic  reproductive  ratio  is  then  7 1,,  —  12/) A'. 


(2.33) 
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2.5  Testing  the  accuracy  of  the  estimation  method 


Before  we  use  our  estimation  on  the  PDV  data,  it  is  important  to  have  a  sense  of 
the  method’s  accuracy  (how  close  are  the  estimates  to  actual,  true  value),  precision 
(how  close  are  the  estimates  from  one  another)  and  potential  bias  (how  far  is  the 
mean  of  the  estimates  from  the  true  value).  I  have  tested  the  method  on  simulated 
epidemic  trajectories  for  a  fixed  infectious  period  equal  to  12  days  over  the  range 
of  initial  susceptible  population  size  So  and  TZq,  So  =  (100;  1,000;  10.000;  100,000), 
and  TZ()  =  ( 1.5, 3,  G,  12).  Figure  2-4  shows  estimates  of  TZo  for  all  the  trajectories  that 
resulted  in  at  least  one  new  infection  (when  there  arc  no  new  infections,  £(p)  -  1  for 
all  values  of  p  so  we  can't  estimate  p  this  way). 

Since  model  (2.24)  is  a  stochastic  one,  there  is  some  probability  not  to  have  an 
outbreak  even  when  one  is  expected.  As  a  result,  when  71q  >  1.  the  distribution  of 
final  sizes  is  bimodal  (as  seen  in  Figure  2-3),  consisting  of  non-outbreaks  and  out¬ 
breaks.  Since  the  realizations  of  the  model  that  are  non-outbreaks  would  not  observed 
in  th<'  available  data  set,  I'm  setting  a  threshold  of  what  constitutes  an  outbreak  so 
as  to  unambiguously  distinguish  between  an  outbreak  and  a  non-outbreak.  I  define 
this  threshold  to  be  20%  of  the  expected  deterministic  final  size  1  —  exp(-/io),  and 
indicate  it  with  a  dotted  line  in  Figure  2-3.  Figure  2-5  summarizes  the  estimates 
from  the  set  of  trajectories  shown  in  Figure  2-4  that  exceed  this  threshold.  The  level 
of  bias  is  comparable  in  both  figures,  but  the  discarding  of  non-outbreak  trajectories 
reduces  the  number  of  outliers,  particularly  for  small  values  of  TZq. 

Each  box  plot  in  Figure  2-5  represents  a  summary  of  estimates  of  IZq  (7 Zo) 
from  1.000  simulated  epidemic  trajectories  with  known  parameters.  I  sealed  the 
value  p  in  each  series  of  simulations  with  respect  to  the  initial  population  size  So. 
p  =  7lo/(12So).  so  that  Ho  is  constant  in  each  graph.  (The  TZU  used  in  simulations 
is  shown  with  horizontal  black  line.)  To  eliminate  the  contributions  of  flic  uncertain¬ 
ties  of  other  parameters,  I  assume  all  other  parameters  (except  p)  are  known  when 
estimating  p  from  the  epidemic  curves. 

The  accuracy  of  the  pseudo-maximum  likelihood  method  described  in  Section  2.4. 1 
depends  on  the  initial  size  of  the  population.  The  method  works  poorly  for  small 
population  sizes  (under  100  individuals),  especially  for  large  values  of  Ho.  Figure  2-5 
also  shows  a  negative  bias;  I  consistently  underestimate  the  true  value  of  7Z{).  which 
is  again  most  prominent  for  small  population  sizes,  and  large  Hq.  For  example,  for 
a  population  of  1000  individuals,  in  78%  cases  of  1  <  Hq  <  1.5,  the  real  value  of  Hq 
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Figure  2-4:  Accuracy  and  precision  of  the  pseudo- maximum- likelihood  method  for 
estimating  7ll}  from  epidemic  curves,  for  various  values  of  IZu  and  and  fixed  in¬ 
fectious  period  (12  days).  Each  box  plot  represents  a  summary  of  estimates  from 
1,000  simulations  using  the  1Z0  value  indicated  by  the  black  horizontal  line.  The  box 
represents  the  inter-quartilc  range  of  the  estimates,  and  the  red  line  is  the  median. 
The  whiskers  are  lines  extending  from  each  end  of  the  box  to  show  the  extent  of  the 
rest  of  the  data;  the  maximum  length  of  the  whiskers  is  1.5  limes  the  inter-quartilc 
range.  Data  that  fall  beyond  the  ends  of  whiskers  are  are  shown  with  red  plus  signs 
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Figure  2-5:  Performance  of  the  pseudo-m aximu m- 1  ikel i  1  lood  method  for  estimating 
from  epidemic  curves,  for  various  values  of  TZ{]  and  Sq.  and  fixed  infectious  period 
{12  days).  Each  box  plot  represents  a  summary  of  estimates  from  1,000  simulations 
using  tlie  71^  value  indicated  by  the  black  horizontal  line.  Trajectories  that  do  not 
exceed  the  threshold  mentioned  in  the  text  are  discarded.  Box  plots  are  constructed 
as  in  Figure  2-4. 
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Figure  2-0:  Each  box  plot  represents  a  summary  of  estimates  from  1. 000  simulations 
using  the  parameters  estimated  from  data,  and  is  constructed  as  in  Fig  2-5.  Box  plots 
on  the  left  show  final  sizes  of  the  simulated  trajectories  used  as  input  epidemic  curves 
to  estimate  1Z{)  values  shown  in  Figure  2-5,  Resulting  estimates  of  the  final  sizes  are 
shown  with  box  plots  in  the  right  column. 
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Table  2.2:  Precision  of  the  pseudo-maximum  likelihood  method.  For  each  range 
below,  simulations  were  made  for  10.000  equally  spaced  values  of  TZd  and  for  s0  = 
1000,  to  =  1.  m  -  0.0.  Trajectories  were  then  used  to  estimate  7?d  The  entries  in  the 
rowrs  of  the  table  fractions  of  simulations  in  each  range  that  resulted  with  a  part  icular 
Hi)  est  imate. 


n  o 

1  -  1.5 

1.5  -2 

2  -  2.5 

2.5  - 

TZq 

3  3-3.5 

3.5-4 

4-4.5 

4.5  -  5 

<i 

1 

0 

0 

0 

0 

0 

0 

0 

1  -  1.5 

0.78 

0.22 

0 

0 

0 

0 

0 

0 

1.5-2 

0.01 

0.79 

0.20 

0 

0 

0 

0 

0 

2  -  2.5 

0 

0.02 

0.74 

0.25 

0 

0 

0 

0 

2.5  -  3 

0 

0 

0.03 

0.65 

0.30 

0.01 

0 

0 

3  -  3.5 

0 

0 

0 

0.05 

0.57 

0.34 

0.04 

0 

>  3.5 

0 

0 

0 

0 

0.02 

0.17 

0.26 

0.55 

comes  from  the  same  range,  and  in  the  other  22%  cases  it  comes  the  range  of  larger 
values.  1.5  <  TZq  <  2  (Table  2.2).  The  negative  bias  is  larger  for  3  <  Tty)  <  2.5, 
when  in  only  57%  of  the  cases  the  true  7£0  lies  in  the  same  range,  34%  of  the  cases 
3.5  <  Hq  <  4.  and  for  4%-  of  the  cases  4  <  TZq  <  4.5.  However,  the  method  does 
extremely  well  for  large  population  sizes  {>  100,000).  for  all  values  of  7Zq  (Figure  2-5. 

One  source  of  bias  is  in  the  round-off  error  in  the  reconstruction  of  the  incidence 
series,  equation  (2.30).  Rounding  off  x  to  the  nearest  integer  towards  minus  infinity 
inevitably  underestimates  the  final  size  of  the  epidemic,  thereby  underestimating  IZq. 
This  bias  is  strongest  for  small  populations  as  the  contribution  of  round-off  error  is 
relatively  larger.  Figure  2-6  compares  the  final  sizes  of  the  epidemic  trajectories  used 
to  estimate  TZ,}  in  Figure  2-5,  and  the  resulting  estimated  final  sizes. 

The  population  sizes  of  almost  all  haul-out  regions  in  the  data  set  that  I  am  using 
to  estimate  of  PDV  outbreaks  fall  in  the  range  1,000  10.000  (see  Table  2.3  for 

details).  For  these  population  sizes,  we  can  expect  our  method  to  give  reliable,  but 
negatively  biased,  estimates  of  TZq. 


2.6  TZq  estimates  from  1988  and  2002  outbreaks 

1  estimated  7ZU  values  for  1988  and  2002  outbreaks  for  the  following  regions:  North 
Skagerrak.  South  Skagerrak,  Swedish  Kattegat,  Danish  Kattegat,  Lim fjord.  Baltic. 
Dutch  Wadden  Sea  (abbreviated  as  WS  NL).  Nieder-Sachsen  Wadden  Sea  (WS  NS). 
Schleswig- Holstein  and  Danish  Wadden  Sea  (WS  SH&DK),  The  Wash,  Tay.  and 


49 


SWEDEN 


Figure  2-7:  Map  of  North  Europe  with  locations  of  PDV  outbreaks  mentioned  in  the 
text. 


Table  2,3:  Population  sizes  of  harbor  seal  populations  before  1988  and  2002  PDV 
outbreaks  for  the  regions  whose  epidemic  curves  were  used  in  estimating  77().  I)(y~) 
designates  the  total  number  of  seals  that  have  died  in  each  location  during  the 
outbreaks  (Harding  e.t  til.,  in  preparation).  Degree  of  spatial  sub-structure  is  given 
by  the  number  of  haul-out  units  that  each  region  consists  of. 


Region 

1988  outbreak 

N  D{o c) 

2002  outbreak 

/V  D(oo) 

haul-out  units 

N  Skagerrak 

2.623 

1.183 

7,466 

4.932 

4 

S  Skagerrak 

n/a 

n/a 

2,702 

1,485 

4 

SW  Kattegat 

2.884 

1 ,823 

4,518 

2,086 

3 

DK  Kattegat 

5.054 

3.266 

0.131 

1.484 

3 

Lim  fjord 

1.474 

614 

1.740 

333 

2 

Balt  ic 

439 

218 

802 

127 

4 

WS  NL 

1 ,800 

914 

7,002 

2,851 

6 

WS  NS 

4.602 

2,634 

10,042 

4,690 

5 

WS  SH & DK 

9,937 

5.774 

18,220 

8.747 

8 

The  Wash 

6.640 

3,535 

6.958 

1,439 

2 

Tay 

700 

31 

1,171 

363 

4 

Moray  Firth 

1.598 

200 

1.198 

86 

5 
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Moray  Firth  (see  map  in  Figure  2-7  for  reference). 

The  population  of  seals  was  unexposed  to  PDV  before  the  1988  outbreak,  as 
indicated  by  serological  studies  (Osterhaus  et  ai.  1988, 1989b;  Thompson  v.t  al.,  1992). 
I  therefore  assume  that  the  entire  population  was  susceptible  at  the  beginning  of  the 
1988  outbreak  (t.  e.,  So  =  N.iq  —  l./’o  =  0),  Exposed  individuals  acquire  life-long 
immunity  to  PDV,  so  a  certain  fraction  of  the  2002  population  consists  of  survivors 
of  the  1988  epizootic.  To  calculate  the  fraction  of  the  2002  population  immune  to 
PDV.  1  assume  that  all  of  the  individuals  were  exposed  to  the  virus  in  1988  (final  size 
=  1),  so  all  of  the  survivors  from  the  1988  epizootic  are  immune.  This  assumption 
is  close  to  the  truth  for  Kattegat  locations,  such  as  Anholt.  where  over  95  %  of  the 
females  were  infected,  as  estimated  from  the  abortion  rates  and  pup  survival  during 
the  epizootic  year  (Heide- Jorgensen  k  Harkbnen,  1992).  The  population  grows  at  the 
rate  A  and  all  of  the  newborn  individuals  are  susceptible  (t.  e..  there  is  no  inherited 
immunity).  Assuming  the  95%  survival  rate  (based  on  data  on  survival  of  adult  seals 
from  Harkbnen  k  Heide- Jorgensen  (1990);  Heide-, Jorgenson  r.t  nl.  (1992))  during  the 
14  years  between  two  outbreaks,  yields  the  fraction  of  the  population  immune  in  2002 
equal  to  0.95 1  x j\u.  Average  population  growth  rates  (A)  for  different  regions  and 
estimates  of  fraction  of  the  population  immune  to  PDV  tire  given  in  Table  2.4. 

In  order  to  reconstruct  the  incidence  series  in  step  (2.30).  we  need  to  know  the 
probability  that  an  infected  individual  dies  from  the  disease  (tn).  This  probability  is 
equal  to  the  ratio  of  the  number  of  individuals  that  have  died  in  the  outbreak.  /?(oc) 
(see  Table  2.3  for  actual  numbers),  and  the  total  number  of  individuals  that  were 
infected.  In  the  absence  of  serological  data,  I  calculated  the  estimates  of  epidemio¬ 
logical  parameters  for  a  range  of  final  sizes  of  the  epidemic:  0.5,  0.0,  0.7,  0.8.  0.9  and 
1.  and  show  the  results  in  Figures  2-8  2-10. 

Final  sizes  below  0.5  are  not  possible  for  most  locations  for  the  observed  values  of 
D(dc).  For  some  localities,  even  some  final  sizes  larger  than  0.5  cannot  be  observed 
for  tin'  number  of  seals  that  have  died  in  those  locations.  In  that  ease,  graphs  for 
these  locations  are  left  blank  for  those  particular  values  of  final  size.  For  example, 
given  the  number  of  seals  that  have  died  in  the  1988  outbreak  in  S\Y  Kattegat,  the 
final  size  of  the  epizootic  for  that  locality  must  have  been  at  least  0.7.  so  there  are  no 
box-plots  for  final  sizes  0.5  and  0.0  in  the  graph  for  SW  Kattegat  in  Figure  2-8. 

To  determine  the  precision  of  the  estimates,  I  simulated  1,000  epidemic  curves 
using  the  estimated  values  of  parameters  for  each  location.  I  discarded  the  non- 
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Table  2,4:  The  fraction  of  the  2002  population  immune  to  PDYr  was  calculated  using 
population  growth  rates  A  and  assuming  95%  survival  rate,  and  the  final  size  of  the 
1988  outbreak  equal  to  1,  Assuming  that  50%  of  susceptible  seals  were  exposed  to 
PDV  in  1988  (final  size  =  0.5)  halves  the  fraction  immune  in  2002. 


Region 

A 

fraction  immune  in  2002 

N.  Skagerrak 

1.14 

0.08 

S.  Skagerrak 

1.13 

0.09 

SW  Katt 

1.13 

0.09 

DKKatt 

1.06 

0.22 

Limfjord 

1.07 

0.19 

Baltic 

1.05 

0.25 

WS  NL 

1.17 

0.05 

WS  NS 

1.11 

0.11 

WS SHDK 

1  13 

0.09 

The  Wash 

1.00 

0.22 

Tay 

1 .04 

0,28 

Moray  Firtli 

0.99 

0.5G 

outbreak  trajectories  (ones  that  did  not  reach  20%  of  the'  approximated  final  size 

1  —  exp(— 7?t|)).  and  estimated  the  values  of  TZq  from  the  remaining  trajectories. 
Figures  2-8  and  2-10  show  the  range  of  estimates  in  box-plots.  Estimates  of  7?,> 
values  fall  in  the  range  1.4  3.15  for  the  1988  outbreak,  and  0.9  3.70  for  the  2002 

outbreak  (over  all  final  sizes). 

2.6.1  Variance  of  7^o  among  locations 

Figures  2-8  -  2-10  illustrate  the  variability  of  estimates  within  a  location.  Estimates 
also  vary  among  locations.  In  this  section  I  address  the  question  of  whether  the  differ¬ 
ence  in  TZq  values  among  locations  is  due  to  estimation  error  or  biological  significance. 

For  the  1988  outbreak,  tin* 1  variance  of  TZq  estimates  assuming  the  final  size  equal 
to  one  is  erf  —  0.103.  Tin1  variance  is  smaller  for  /  =  0.7,  <7q7  =  0.014.  Can  this 
variance  in  estimates  Ire  observed  by  chance  alone  (the  null  hypothesis  //0),  or  are 
the  differences  in  7ZU  among  locations  statistically  significant?  To  find  out,  I  did  a 
randomization  test  by  permuting  observations  among  locations  (c.  <j.  Caswell.  2001). 
For  the  1988  outbreak,  the  data  consists  of  11  epidemic  curves,  and  I  treat  each 
recovered  seal  with  its  respective  relative  day  of  recovery  (day  since  the  first  seal  was 
found  dead  in  that  location)  as  a  data  point.  Under  //u,  the  combination  of  time- 
series  of  dead  seals  that  consist  the  observed  epidemic  curves  for  different  locations 
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Figure  2-8:  Estimates  of  7Zq  from  data  (black  dots)  assuming  different  filial  sixths  of  t  lie 
epidemic  for  various  locations  in  Europe  for  1988  outbreak.  Each  box  plot  represents 
a  summary  of  estimates  from  1,000  simulations  using  1  lit'  parameters  estimated  from 
data,  and  is  constructed  as  in  Fig  2-5. 
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Figure  2-U:  Box  plots  summarize  estimates  of  the  final  size  of  the  1,000  trajectories 
simulated  with  the  parameters  corresponding  to  Figure  2-8  and  final  size  1, 
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Figure  2-10:  Estimates  of  TZq  from  data  (black  dots)  assuming  different  tinal  sizes 
of  the  epidemic  for  various  locations  in  Europe  for  2002  epizootic.  A  fraction  of 
the  population  that  is  assumed  immune  in  final  size  “1  im.”  is  shown  in  Table  2.4. 
Each  box  plot  represents  a  summary  of  estimates  from  1,000  simulations  using  the 
parameters  estimated  from  data,  and  is  constructed  as  in  Fig  2-5. 
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Figure  2-11:  Randomization  test  for  the  1988  outbreak  assuming  two  different  values 
of  final  size.  /:  A)  f  =  1:  B)  /  =  0.7.  Dotted  lines  represent  the  variance  in  7lu  from 
original  data,  of  =  0.103  and  <t£7  —  0.014. 

is  just  as  likely  as  any  other  series,  and  can  be  observed  by  chance  alone.  To  obtain 
the  distribution  of  the  variance  of  it ()  under  //,).  1  permuted  the  seals  to  create  new 
epidemic  curves  for  each  location,  by  randomly  drawing  /t,(oc)  individual  dead  seals 
from  the  original  data  set  without  replacement,  (/tj(oc)  is  the  total  number  of  seals 
that  have  died  in  location  /  given  in  Table  2.3.)  This  is  repeated  10.000  times  and 
for  each  of  the  10,000  permutations,  I  estimate  for  all  locations,  calculate  the 
variance  of  the  estimates,  and  obtain  a  distribution  of  variance  shown  in  Figure  2-11. 
For  /  =  1,  the  probability  of  observing  the  variance  in  7^o  that  is  larger  than  one 
observed  in  the  original  data  by  chance  alone  is  about  1  in  10,000.  For  /  =  0.7,  the 
probability  of  observing  larger  variance  than  in  the  original  data  is  less  than  0.04. 

The  randomization  test  for  the  2002  outbreak,  shows  that  the  variance  in  values 
falls  outside  of  the  distribution  of  variance  under  //().  The  variance  in  7£()  estimates 
assuming  final  sizes  /  =  1  and  /  =  0.7  are  af  =  (J.35  and  —  0.07.  respectively. 
Figure  2-12  shows  that  the  probability  of  observing  the  variance  greater  than  or  equal 
to  one  observed  in  data  by  chance  alone  is  less  than  1  in  10,000.  Therefore,  we  can 
conclude  that  the  observed  difference  in  estimates  among  locations  is  statistically 
significant. 

Since  the  variance  in  7vn  is  significant,  the  observed  difference  is  a  result  of  a 
biological  or  epidemiological  process,  or  habitat  differences.  To  identify  the  patterns 
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Figure  2-12:  Randomization  test  for  the  2002  outbreak  assuming  two  different  values 
of  final  size,  /:  A)  /  =  1;  B)  /  =  0.7.  Dotted  lines  represent  the  variance  in  72(i  from 
original  data,  a'f  =  0.35  and  <7q7  =  0.07. 
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Figure  2-13:  Comparison  of  regional  IZq  estimates  between  two  outbreaks. 
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Figure  2-14:  Relationship  of  TZ()  and  the  pre-epizootic  population  size. 
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Figure  2-15:  Relationship  of  TZ{)  and  duration  of  the  epizootic  (in  days). 
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Figure  2-16;  The  proportion  of  the  population  killed  during  the  outbreak  decreases 
with  the  peak  mortality  dateXso.  the  day  of  the  year  when  50 %  of  the  final  disease 
mortality  was  reached. 
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Figure  2-17;  Relationship  of  TZq  and  Xj0. 
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Figure  2-18:  The  influence  of  topography  on  estimates. 


Figure  2-19:  Epidemic  curves  for  2002  outbreak  for  3  haul-out  locations  in  Danish 
Kattegat,  and  pooled  data  for  Danish  Kattegat.  Dotted  lines  represent  first  5%  and 
last  5%  of  the  reported  eases.  The  numbers  in  the  parentheses  give  the  duration  of 
each  outbreak  in  days,  defined  as  the  period  between  the  first  5%  and  the  last  5%  of 
the  reported  cases  (indicated  by  dotted  lines). 
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Figure  2-20:  Relationship  of  TZq  and  the  number  of  haul-out  units  within  each  region 
for  1988  and  2002  outbreaks. 
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Figure  2-2 L:  Relationship  of  effective  reproductive  number,  TZr,  and  population 
growth  rate,  and  fraction  immune  for  2002  outbreak. 
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in  7l0  differences  and  what  influences  them,  I  examine  the  relationship  of  TZo  to  initial 
susceptible  population  size,  duration  of  the  epizootic,  peak  mortality  date,  degree  of 
spatial  structure,  and  the  level  of  immunity. 

For  density-dependent  transmission,  we  expect  7Zo  to  be  proportional  to  popula¬ 
tion  size,  all  other  parameters  being  equal.  Since  the  population  level  on  the  continent 
in  2002  was  about  double  that  of  1088.  wc  would  expect  to  observe  higher  7llt  values 
on  1  he  continent  in  2002  than  in  1988.  assuming  that  all  else  were  equal.  However.  7Z{, 
estimates  shown  in  Figure  2-13  show  no  consistent  differences  between  two  outbreaks, 
and  TZu  does  not  depend  on  the  population  size  (Figure  2-14). 

Tiie  relationship  of  TZq  to  the  duration  of  the  epizootic  is  shown  in  Figure  2- 
15.  Differences  in  sampling  and  reporting  effort  are  highest  at  the  beginning  and 
the  end  of  the  epizootic,  so  I  represent  the  duration  of  the  epizootic  (AT)  with  the 
period  between  the  first  5%  and  last  5%  of  the  dead  seals  reported  for  each  location 
(Harkonen  et  al.  2006).  The  higher  7£(J  the  faster  disease  will  spread  and  the  shorter 
its  duration. 

Phoeine  distemper  peaked  in  different  times  of  the  year  in  different  regions.  1 
describe  this  different  timing  by  7},o.  or  the  Julian  day  of  the  year  when  half  of  the 
final  mortality  in  a  particular  region  was  reached.  This  parameter  is  more  reliable 
than  the  date  of  the  first  finding  (the  day  when  the  first  dead  seal  was  reported  in 
each  region)  as  there  is  high  uncertainty  in  the  reporting  at  the  beginning  of  the 
outbreak.  In  general,  the  fraction  killed  by  PDV  decreases  with  T™  and  is  smallest 
for  locations  where  the  outbreak  started  at  the  latest  (Figure  2-16).  For  tin1  1988 
outbreak,  the  value  of  TZ^  drops  in  the  regions  where  the  outbreak  peaked  late  in 
the  year  (Figure  2-17).  suggesting  that  there  might  be  a  seasonal  mechanism  that 
affects  the  transmission  of  the  virus.  For  2002  outbreak,  this  relationship  is  obscured 
by  tht'  Limfjord  locality.  The  Limfjord  population  suffered  an  additional  source  of 
mortality  in  2002,  linked  to  malnutrition  (Karin  Harding,  personal  communication), 
that  can  increase  the  observed  7^o-  It  has  been  difficult  to  discriminate  disease- related 
mortality  from  other  causes  of  death,  so  the  contribution  of  the  non- disease- related 
mortality  to  the  value  of  the  Limfjord  7?n  remains  unknown. 

Timing  of  infection,  together  with  the  topography  of  the  location,  can  influence 
the  dynamics  of  an  outbreak.  Figure  2-18  groups  the  locations  according  to  their 
sea-bottom  substrate,  into  rocky  or  sandy  regions.  In  the  rocky  regions,  7 drops 
with  the  peak  mortality  date,  whereas  in  the  sandy  regions  this  relationship  is  absent. 
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Spatial  structuring  of  the  epidemic  data  can  also  influence  the  estimate  of  7£0, 
PDV  sets  up  a  sequence  of  local  outbreaks,  that,  once  initiated,  proceed  indepen¬ 
dently.  Pooling  data  from  several  localities  from  a  larger  region  into  a  single  epidemic 
curve,  prolongs  the  duration  of  the  epizootic  in  that  region  and  decreases  the  slope 
of  the  epidemic  curve,  thereby  decreasing  1 I  illustrate  this  with  Danish  Kattegat 
data  from  the  2002  outbreak  in  Figure  2-19,  for  which  the  individual  epidemic  curves 
from  different  haul-out  units  are  available  (Harkonen  et  al..  2000).  The  larger  the 
number  of  the  haul-out  units  with  different  starting  dates  of  the  outbreak  included 
in  the  epidemic  curve,  the  smaller  7Z0.  This  expected  relationship  is  also  observed 
for  the  1988  outbreak  (Figure  2-20).  and  7£()  drops  as  the  spatial  structuring,  i. 
the  number  of  haul-out  units,  in  a  particular  region  increases.  The  pattern  disap¬ 
pears  in  2002.  suggesting  that  another  process,  such  as  immunity,  absent  in  1988.  is 
influencing  7£»  in  2002. 

A  certain  fraction  of  the  2002  population  are  the  survivors  from  the  1988  epi¬ 
zootic.  It  is  therefore  natural  to  assume  that  immunity  could  play  an  important  role 
in  the  2002  outbreak,  and  change  the  patterns  observed  in  the  1988  outbreak.  To 
account  for  immunity.  I  have  estimated  the  fraction  of  each  population  that  can  be 
immune  to  PDV  in  2002  (Table  2.4).  and  have  estimated  new  7ZU  values  by  decreasing 
initial  susceptible  population  accordingly.  Correcting  for  immunity  did  not  change 
the  estimates  of  TZ{)  significantly  (Figure  2-10).  because  in  general  the  fraction  of  the 
population  that  was  likely  to  be  immune  in  2002  is  small. 


2.7  Discussion 

Estimating  epidemiological  parameters  from  data  is  a  challenging  task.  Information 
on  infectious  diseases  is  frequently  lacking  in  detail,  since  it  is  hard  to  observe  the 
exact  times  when  individuals  become  infected  and  by  whom.  Data  on  infectious 
diseases  in  wildlife  are  even  more  incomplete  than  human  epidemiological  data  or 
data  on  diseases  in  domestic  animals,  In  the  case  of  phocine  distemper  the  only 
available  data  is  the  number  of  seals  that  have  died  from  the  disease, 

I  present  a  novel  method  for  estimating  that  can  be  useful  for  epidemics  where 
the  number  of  dead  is  the  only  information  we  observe.  Model  (2.24)  assumes  a 
closed  epidemic,  so  the  method  is  not  applicable  for  diseases  with  long  infectious 
periods,  such  as  HIV/AIDS,  where  the  susceptible  population  is  being  replenished 
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by  birth  during  the  outbreak.  Even  though  the  method  is  ad  hoc  I  tried  to  base  it 
on  the  maximum-likelihood  framework.  Simulation  results  show  that  the  estimates 
are  reliable  for  population  sizes  in  the  range  covered  by  the  data.  The  method  is 
negatively  biased,  but  for  the  observed  population  sizes  this  bias  is  small.  One  source 
of  bias  is  numerical,  due  to  the  round-off  error  in  the  reconstruction  of  the  time-series 
for  incidence.  Additional  sources  of  bias  can  come  from  assumptions  about  infectious 
period  and  from  observational  errors.  We  assume  a  3-day  latent  period  followed  by 
a  12-day  infectious  period.  In  reality,  infectious  period  can  vary  between  5  and  10 
days  (Harder  fit  al. .  1990).  which  can  affect  overall  transmission  during  the  period  of 
infection  and  thereby  influence  the  estimated  number  of  secondary  cases.  Another 
source  of  error  is  "built  into”  our  series  of  dead  seals,  which  we  use  to  reconstruct  the 
time  series  of  incidence!  and  infectious  case  counts.  Data  on  seals  that  have  died  from 
PDV  is  collected  by  counting  carcasses  that  stranded  ashore,  which  depends  on  the 
weather  conditions  and  sampling  and  reporting  effort,  so  there  is  inevitably  an  error 
associated  with  the  observation  process. 

Averaging  data  over  large-scale  spatial  structure  can  lead  to  further  underesti¬ 
mation  of  71q.  Phocine  distemper  epizootic  sets  up  a  sequence  of  local  outbreaks 
that,  once  initiated,  proceed  independently.  The  spatial  unit  at  which  an  outbreak 
occurs  is  a  single  haul-out  location.  Treating  several  haul-out  units  as  one,  or  pool¬ 
ing  data  from  several  haul-out  units  into  one,  will  result  in  a  smaller  estimate  of 
than  for  a  single  haul-out.  Looking  at  the  outbreak  at  a  larger  spatial  scale  than  it 
occurs,  artificially  prolongs  the  duration  of  the  epizootic  thereby  reducing  R{)  (see 
Figure  2-15). 

One  potential  way  to  improve  the  method  and  reduce  the  bias  would  be  to  use  EM 
algorithm  to  find  the  maximum-likelihood  estimates.  The  EM  algorithm,  first  named 
by  Dempster  et  al.  (1977)  is  a  computation  that  iterates  between  an  "Expcctation- 
step"  and  a  "Maximization-step.”  It  is  a  broadly  applicable  algorithm  for  estimating 
parameters  from  incomplete  data  sets  by  maximum  likelihood  methods.  In  epidemi¬ 
ology.  the  EM  algorithm  is  not  only  used  for  the  estimation  of  parameters,  but  also 
for  reconstruction  of  the  entire  time  series  of  unobservable  classes.  For  example,  the 
most  reliable  data  on  HIV/ AIDS  epidemic  are  often  incidences  of  eases  diagnosed 
with  AIDS.  It  is  much  harder  to  observe  the  times  of  infection  with  HIV.  Becker 
(1997)  used  the  EM  algorithm  to  back-project  HIV  infection  curve  from  the  AIDS 
incidence  data.  Andersson  Britton  (2000)  show  how  the  algorithm  can  be  used 
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to  estimate  the  probability  of  infection  in  chain-binomial  models.  In  the  case  of  the 
distemper  outbreak  in  seals,  the  only  time-series  data  at  hand  is  actually  the  series 
of  estimates  of  the  number  of  seals  that  have  died  each  day,  so  EM  approach  is  not 
applicable. 

Estimates  of  720  values  for  phocine  distemper  outbreaks  fall  in  the  range  1.4 
3.15  for  the  1988  outbreak,  and  0.9-  3. 70  for  the  2002  outbreak,  over  all  final  sizes. 
Since  the  method  described  in  Section  2.3  is  negatively  biased,  the  true  value  of  72« 
for  phocine  distemper  outbreaks  is  likely  to  be  higher.  Estimates  for  1988  outbreak 
reported  in  this  chapter  agree  with  72, }  for  the  1988  epizootic  from  the  literature  72 
was  estimated  to  be  2.8  by  Swinton  et  al  (1998)  and  between  2.1  and  3  by  de  Koeijer 
et  al.  (1998).  There  are  no  published  estimates  of  720  for  the  2002  outbreak. 

Phocine  distemper  virus  has  a  small  72o  compared  to  some  of  its  relatives  be¬ 
longing  to  the  genus  Morbilliviridae.  Measles,  one  of  the  most  infectious  childhood 
diseases  and  the  most  famous  member  of  the  Morbilli virus  group,  has  72,,  that  ranges 
between  10  and  20.  depending  on  location  (Anderson,  1990;  Keeling  k  Grenfell,  2000; 
Bjarnstad  et.  al.  2002),  Highly  transmissible  infections  include  some  other  childhood 
diseases  like  chickcupox  (7  <  72o  <  12).  mumps  (11  <  72,)  <  14),  and  also  sexually 
transmitted  diseases  such  as  HIV/ AIDS  epidemic  with  72o  around  10  in  sub-Saharan 
Africa  (Anderson,  1990).  However,  even  diseases  with  small  72,,  and  transmissibility 
compared  to  measles  can  cause  pandemics  and  outbreaks  of  severe  morbidity  and 
mortality;  estimates  of  72o  for  1918  pandemic  influenza  lie  between  2  and  3  (Mills 
et  al. ,  2004).  whereas  seasonal  influenza  epidemics  have  72(,  around  1.35  (Viboud 
et  al..  2006),  and  SARS  outbreak  of  2002  had  72t)  value  around  3  (Lipsitch  et  al.. 
2003)  In  naive  populations,  diseases  with  relatively  low  transmissibility  and  72a  can 
cause  severe  morbidity  and  mortality,  which  appears  also  to  have  been  the*  case  with 
phocine  distemper  outbreaks  in  harbor  seal  populations. 

Values  of  72o  estimated  for  PDV  can  also  give  us  insight  on  the  transmission 
process  going  on  between  seals.  If  transmission  is  density-dependent,  as  has  been 
assumed  in  the  model  (2.24),  72,,  should  increase  with  initial  population  size,  ,S’0,  all 
else  being  equal.  In  the  case  with  phocine  distemper,  there  is  no  clear  relationships 
between  Sq  and  72a  (Figure  2-14).  Since  Figure  2-14  compares  different  locations, 
the  ‘all  else  being  equal’  assumption  does  not  hold,  so  it’s  hard  to  determine  how 
much  other  habitat  differences  contribute  to  relationship  of  So  and  72o.  Nevertheless, 
estimates  of  720  are  not  very  different  (they  arc  all  the  same  order  of  magnitude) 


even  though  populations  sizes  varied  between  400  and  20,000,  On  the  other  hand, 
for  frequency-dependent  transmission  TZq  in  equation  (2,6)  is  constant  irrespective  of 
So.  This  suggests  that  in  reality  the  mixing  of  seals  is  somewhere  in  the  continuum 
between  those  to  extremes. 

Processes  that  affect  the  mixing  of  seals  are  likely  to  affect  the  value  of  Ho  as  well. 
Figure  2-17  suggest  that  Ho  is  smaller  for  the  locations  where  the  disease  appeared 
late  in  the  year,  especially  in  the  1988  outbreak.  Seasonal  processes  that  result  in 
lower  mixing  rates  late  in  tin1  year  can  help  explain  this  pattern.  One  of  those  seasonal 
processes  is  the  haul-out  behavior  of  seals.  Seals  give  birth,  rear  their  offspring  and 
molt  on  land,  leading  to  their  haul-out  behavior  to  peak  in  the  summer,  when  up 
to  67%  of  the  colony  can  be  found  on  land  (Heide-, Jorgensen  k  Harkonen.  1992). 
Since  PDV  is  an  airborne  virus,  it  is  believed  its  spreads  by  inhalation  while  seals 
are  hauled-out  on  land.  Seasonal  haul-out  behavior  can  influence  the  mixing  and  the 
contact  processes  in  seals  leading  to  reduced  transmission  of  infection  in  the  fall  and 
winter,  and  to  small  Hq  values. 

The  haul-out  behavior  differs  from  region  to  region,  and  in  Kattegat  and  Skagcrrak 
it  depends  on  the  sea-bottom  substrate  of  the  locations.  This  also  influences  the 
dynamics  of  the  outbreak,  so  in  the  rocky  regions  Ho  decreases  with  the  peak  mortality 
date,  whereas  in  the  sandy  regions  this  is  not  so.  In  the  rocky  regions,  food  is  abundant 
and  seals  don’t  have  to  travel  far  in  search  for  food,  so  they  spend  more*  time  on  land  . 
In  the  sandy  regions,  the  food  is  scarce  so  seals  spend  less  time  on  land,  and  more 
time  in  search  for  food.  Throughout  the  year  the  fraction  of  the  population  hauled- 
out  on  land  is  lower  in  the  sandy  regions  than  in  the  rocky  regions,  so  it  will  be  more 
influenced  by  the  stochasti cities  of  the  epidemic  processes  leading  to  more  variable 
H0  estimates.  The  effects  of  the  seasonal  lmul-out  behavior  and  the  timing  of  the 
outbreak  on  the  dynamics  of  the  phoeine  distemper  epizootics  are  studied  in  further 
detail  in  Chapter  3. 
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Chapter  3 


Seasonal  haul-out  behavior  and  the 
dynamics  of  the  phocine  distemper 
virus 

3.1  Introduction 

Phocine  distemper  virus  was  first  described  in  1988  after  it  killed  over  23.000  harbor 
seals  (Phoca  minima)  in  Northern  Europe  (Osterhaus.  1988;  Osterhaus  A;  Vedder, 
1988;  Osterhaus  et  al. ,  1989a;  Dietz  et  al.,  1989;  Hekle-Jorgensen  et.  al.,  1992).  The 
disease  appeared  at  the  Danish  island  of  Anholt  in  April,  and  in  the  following  months 
spread  throughout  Europe,  reaching  Scotland  in  the  fall.  Over  the  next  14  years  the 
harbor  seals  populations  recovered  to  double  the  1988  pre-epizootic  population  size 
on  the  Continent,  and  reaching  levels  comparable  to  1988  in  the  UK. 

An  interesting  feature  of  PDV  outbreaks  is  that  the  proportion  of  seals  that  died 
during  tin*  outbreak  varied  among  regions  (Figure  3-1).  The  Danish.  Swedish  and 
Norwegian  populations  experienced  much  higher  mortalities  (50  —  00%)  than  the 
populations  in  England,  Scotland  and  Ireland  (10—20%)  (Dietz  et  al.,  1989;  de  Koeijcr 
et  al..  1998;  Harding  et  al. ,  2002;  Harkonen  et  al..  2000).  The  mortality  was  lowest  in 
the  regions  where  PDV  appeared  last.  Estimates  of  Ra  also  differ  significantly  among 
locations  (Chapter  2).  Rq  was  estimated  to  be  lower  in  locations  where  the  disease 
appeared  late  in  the  year,  suggesting  a  seasonal  mechanism  might  be  influencing  the 
transmission  of  the  virus.  One  potential  mechanism  is  the  roughly  annual  cycle  in 
the  haul-out  behavior  of  harbor  seals. 
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Figure  3-1:  Map  of  mortalities  of  1988  and  2002  PDV  outbreaks. 

Seals  give  birth,  rear  their  offspring  and  molt  on  land.  As  a  result,  the  fraction 
of  a  population  "hauled-out"  on  land  at  any  time  peaks  fluring  the  pupping  and 
molting  seasons  in  late  spring  and  summer.  At  that  time,  approximately  00%  of  a 
colony  are  on  land  {Thompson  <:f  ai.  1997;  Ries  rt  ai,  1998;  Harkoncn  r  f  ai,  1999). 
The  percentage  of  animals  on  land  also  varies  with  age  and  sex  (Thompson,  1989; 
Harkonen  et  ai.  1999.  2002). 

In  addition  to  demography,  the  haul-out  pattern  is  influenced  by  the  sea-bottom 
topography.  In  tin'  regions  of  the  Kattegat  and  Skagerrak  Seas  with  rocky  sea-bottom 
substrate,  food  is  abundant  and  seals  do  not  have  to  travel  far  in  search  for  food,  so 
they  spend  more  time  on  land.  Where'  the  substrate  is  sandy,  food  is  scarce  so  seals 
Spend  less  time  on  land,  and  more  time  in  search  for  food.  As  a  result,  throughout  the 
year  the  fraction  of  the  population  hauled-out  on  land  is  lower  in  the  sandy  regions 
than  in  the  rocky  regions  (Figure  3-2). 

Phocine  distemper  virus  is  ail  airborne  virus,  that  spreads  by  inhalation  (Kennedy. 
1990,  1998)  and  can  only  be  transmitted  between  seals  that  arc  hauled-out  on  laud 
Haul-out  behavior  determines  the  contact  process  between  the  seals,  so  it  will  influ¬ 
ence  the  transmission  of  the  virus  and  make  it  seasonal  whenever  haul-out  behavior 
is  seasonal.  In  this  chapter.  I  account  for  the  haul-out  behavior  by  including  it  in 
t  lit*  model  of  transmission.  Tin1  result  is  a  noil-autonomous  model  for  the  epidemic. 
I  modify  the  pseudo-maximum-likelihood  estimation  procedure  from  Chapter  2  to 
estimate  the  probability  of  infection,  and  derive  an  expression  for  72u. 
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sandy  substrate 


rocky  substrate 


Fi^urt‘  3-2:  Haul-out  patterns  in  Kattegat  (left)  and  Skagcrrak  Seas  (right)  depend 
on  sea  bottom  topography. 

When  the  model  includes  the  haul-out  behavior,  the  importance  of  the  timing  of 
the  virus  introduction  becomes  apparent.  1  use  simulations  to  illustrate  how  seasonal 
behavior  and  the  timing  of  the  beginning  of  infection  influence  the  dynamics  of  the 
outbreak.  Mv  results  show  that  the  mortality  and  tin'  final  size  of  the  outbreak  will 
be  low  if  the  virus  is  introduced  to  the  population  in  the  winter,  when  tht'  population 
numbers  on  land  are  lowest . 


3.2  Model 

The  duration  of  phocine  distemper  outbreaks  is  short  compared  to  the  lifespan  of 
harbor  seals,  so  I  assume  that  the  seal  population  did  not  grow  significantly  during 
any  of  the  two  epizootics  and  model  PDV  outbreaks  in  any  haul-out  location  as 
a  closed  epidemic.  The  seal  population  on  day  /  consists  of  susceptible  seals  (St). 
infectious  seals  (/().  and  a  removed  class  (Rt).  which  accounts  for  both  immune  and 
dead  seals.  The  incidence,  or  the  number  of  seals  that  are  infected  on  day  /.  is  given 
by  the  random  variable  Xt.  The  random  variables  are  denoted  by  upper-case  letters 
here,  and  their  realizations  by  lower  ease.  After  getting  infected,  a  seal  goes  through 
a  latent  period  of  3  days,  after  which  it  becomes  and  remains  infectious  For  12  days, 
when  it  finally  becomes  immune  or  dies.  In  Chapter  2.  I  modeled  the  dynamics  of 


tliis  epidemic  process  as 


i+i 

=  &i  ~  Aj, 

(3.1a) 

HI 

=  it  +  ^t-3  — 

(3  lb) 

'HI 

=  rt  +  Xt-io- 

(3.1c) 

Xt 

~  Binfsi,  1  —  (1  —  p)*e], 

(3. Id) 

where  xf  is  zero  for  i  negative  (Heide- Jorgensen  <k  Harkduen,  1 U02 } , 

System  (3.1)  does  not  account  for  haul-out  behavior,  so  a  contact  between  any  pair 
of  infectious  and  susceptible  seals  can  result  in  a  new  infection,  with  the  probability 
p.  The  probability  that  a  susceptible  docs  not  get  infected  after  contacting  a  given 
infective  is  1  —p.  and  (1  —  p)lt  is  the  probability  of  avoiding  getting  infected  by  any 
of  the  it  infect ives  at  time  /.  The  total  probability  that  a  susceptible  gets  infected  on 
day  l  is  then  1  —  (1  — 

Since  PDV  is  an  airborne  virus.  I  assume  that  it  can  only  be  transmitted  between 
seals  that  are  hauled-out.  On  any  given  flay  f,  every  seal  lias  the  same  probability 
fit  of  being  hauled-out.  Let  and  l_t  be  biliomially  distributed  random  variables 
that  describe  the  number  of  susceptible  and  infectious  seals  hauled-out  on  day  I. 
respectively. 


St  ~  Bm[sf,ht],  (3.2a) 

/,  ~  Bin [it (3.2b) 

The  probability  that  an  infective  seal  on  land  meets  and  infects  a  susceptible  seal 
on  land  during  one  day  is  p.  The  incidence,  A',,  is  now  a  function  of  hauled-out 
susceptibles  and  hauled-out  infectious  seals 

Xt  ~  Bui[st,l-  (l-p)H  (3.3) 

After  the  infectious  period,  seals  either  recover  or  die.  If  the  the  probability  of  death 
is  m.  the  number  of  seals  that  die  on  day  (  (V))  is  biliomially  distributed 

Y,  ~  Bin[xt_i5,m].  (3.4) 

System  (3.1)  (3.3)  is  non-aUtonomous;  the  probability  of  being  hauled-out  on  land  is 
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a  function  of  time. 


3.2.1  7£o  for  constant  haul-out 

Seasonal  behavior  of  seals  affects  the  dynamics  of  the  model.  We  are  particularly 
interested  in  how  the  haul-out  changes  the  threshold  behavior  of  the  mode]  and  its 
7£().  In  nature,  the  seasonal  behavior  of  seals  varies  in  time,  but  for  mathematical 
simplicity  I  will  first  derive  an  expression  for  7?0  for  the  case  where  haul-out  is  constant 
in  time. 

When  the  fraction  of  seals  on  land  is  constant  in  time,  ht  -  h,  seals  have  a  fixed 
probability  of  being  hauled  out  on  land  throughout  the  year.  For  a  given  j>.  the 
number  of  seals  that  become  infected  during  the  outbreak  depends  on  the  value  of 
h  the  smaller  value  of  h.  the  smaller  the  size  of  the  outbreak  (Figure  3-3}.  When 
h  =  1,  the  model  reduces  to  the  model  without  haul-out  behavior  from  Chapter  2. 
given  in  equations  (2.24)  and  (2,25). 

The  basic  reproductive  number  now  depends  on  the  haul-out  behavior  in 
addition  to  transmission  probability.  One  infectious  individual  in  a  population  of  so 
suseept idles,  will  on  average  cause  (pkst))h  new  infections  during  one  day.  as  it  has 
the  probability  of  being  on  land  equal  to  h,  and  once  it  is  on  land  it  will  come  into 
contact  with  hsQ  susceptibles.  An  individual  remains  infectious  for  12  days,  so  the 
expected  number  of  new  infections  it  will  produce  during  its  infectious  lifetime  is 
equal  to 

1Z0  =  12/>Vo-  (3.5) 

For  Tli)  <  1  only  small  outbreaks  outbreaks  occur,  and  the  distribution  of  final 
sizes  is  unimodal.  When  7^()  >  1  both  minor  and  major  outbreaks  are  possible,  and 
the  final  sizes  have  a  bimodal  distribution  (Figures  3-3  and  3-4). 

3.2.2  7Z{)  for  time- varying  haul-out  behavior 

When  haul-out  behavior  is  constant  (or  absent ) ,  the  probability  of  an  outbreak  is 
the  same  regardless  of  what  day  the  virus  is  introduced  to  the  population.  Equal 
percentage  of  the  population  is  on  land  every  day,  so  for  every  day  of  the  virus  in¬ 
troduction.  the  expected  number  of  new  infections  is  the  same.  TZq  does  not  change 
with  the  timing  of  the  virus  introduction.  When  different  percentages  of  the  popula¬ 
tion  are  hauled-out  at  different  times  of  the  year,  the  expected  number  of  seals  that 
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Figure  3-3:  The  final  size  of  an  outbreak  depends  on  the  fraction  of  seals  hauled 
out  on  land,  h.  Parameter  values:  So  =  1000,  p  =  0.0005,  m  =  0.0.  The  TZu  values 
corresponding  to  these  parameter  values  are  indicated  on  the  graph. 
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Figure  3-4:  The  distribution  of  final  sizes  as  the  function  of  the  haul-out  fraction. 
For  each  of  the  values  of  h.  this  graph  shows  final  sizes  of  1000  epidemic  trajectories 
with  equal  parameter  values  {/;  =  0.0005,  —  1000).  The  vertical  line  designates 

the  value  of  haul-out  fraction  for  which  Hq  —  1.  When  7Z[}  >  1  major  outbreaks  are 
possible,  and  the  distribution  of  final  sizes  is  bi  modal. 

become  infected  depends  on  the  day  the  virus  is  introduced  to  the  population.  When 
the  number  of  seals  on  land  is  low,  the  potential  for  a  successful  contact  is  small  and 
there  is  no  outbreak.  As  the  number  of  seals  on  land  increases,  so  does  the  probability 
and  final  size  of  outbreak. 

Figure  3-5  shows  1,000  epidemic  trajectories  for  three  different  days  of  the  in¬ 
troduction  of  infection.  In  Fig  3-5A.  one  infectious  individual  is  introduced  to  an 
otherwise  susceptible  population  on  Julian  day  110.  when  only  about  10%  of  the  pop- 
ulation  is  hauled-out  on  land.  In  this  ease,  a  major  outbreak  occurs  in  only  4  of  tin1 
1.000  trajectories.  If  we  introduce  one  infectious  seal  to  the  population  on  Julian  day 
ICO,  almost  every  simulation  results  in  an  outbreak,  but  the  total  number  of  seals 
infected  during  the  outbreak  (i,  e.,  the  final  size)  varies.  If  the  disease  appears  later 
in  the  year,  Julian  day  200.  when  the  numbers  on  land  are  declining,  the  final  size  of 
the  outbreak  decreases,  as  does  the  final  mortality, 

Figure  3-5  shows  only  3  of  the  possible  305  days  on  which  the  infection  may  begin. 
To  illustrate  how  the  mortality  and  the  final  size  of  the  epidemic  would  change  with 
different  starting  days  of  the  infection.  I  introduce  one  infectious  seals  to  an  otherwise 
susceptible  population  on  each  day  of  the  year.  For  /p,  I  use  two  different  types  of 
haul-out  curves  shown  in  Figure  3-2.  for  sandy  and  rocky  sea-floor  topography.  For 
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Figure  3-5:  Influence  of  the  day  of  the  introduction  of  the  virus  on  the  dynamics  of 
the  outbreak.  Graphs  on  the  top  show  the  haul-out  pattern  for  the  sandy  regions 
that  was  used  in  simulations.  The  vertical  lines  in  the  top  graphs  indicate  the  day 
that  the  first  infectious  seal  appeared  (to).  Each  graph  on  the  bottom  shows  1.000 
epidemic  trajectories  for  given  above.  All  other  parameters  are  equal  in  the  graphs 
S0  =  1500.  p  =  0.001,  m  =  0.0.  (A)  t.0  =  110,  (B)  L0  =  100.  (C)  t0  =  200. 
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day  of  virus  introduction 


Figure’  3-0:  Seasonal  haul-out  behavior  and  the  timing  of  the  virus  introduction 
both  influence  the  mortality  observed  at  the  end  of  the  outbreak.  For  each  day  of 
the  year  I  simulated  10.000  epidemic  trajectories  using  that,  day  as  the  day  of  virus 
introduction.  I  summarize  the  final  mortality  and  final  siy,e  (total  number  infected 
during  an  outbreak)  observed  in  those  10,000  simulation  with  the  mean  (solid  line) 
and  with  the  5th  and  95th  percentile  (gray  envelope).  Sq  =  1500,  ph  —  0.005.  m  =  0.6. 
A)  sandy  regions,  B)  rocky  regions. 
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Figure  3-7:  Seasonal  haul-out  behavior  and  timing  of  the  virus  introduction  influence 
the  mortality  observed  at  the  end  of  the  outbreak,  For  each  day  of  the  year  1  simu¬ 
lated  10,000  epidemic  trajectories  using  that  day  as  the  day  of  virus  introduction,  1 
summarize  the  final  mortality  and  final  size  observed  in  those  10,000  simulation  with 
the  mean  (solid  line)  and  with  the  5th  and  05th  percentile  (gray  envelope).  Su  =  1500, 
pi,  =  0.001.  m  -  0-6.  A)  sandy  regions.  B)  rocky  regions. 
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each  day  of  the  virus  introduction.  I  simulate  10.000  epidemic,  trajectories  300  days 
long.  The  last  point  of  each  trajectory  provides  the  final  mortality.  For  each  starting 
day  of  the  outbreak,  I  get  10,000  values  of  the  final  mortality.  I  summarize'  those 
values  of  final  mortality  for  different  values  of  p  in  Figures  3-0  and  3-7  with  the  mean 
value  (solid  line),  and  the  envelope  defined  by  the  5th  and  95th  percentile. 

In  the  winter  (the  beginning  and  the  end  of  the  year;  approximate  Julian  days 
0-100.  350-365).  there  is  virtually  no  possibility  of  an  outbreak  as  only  a  small  number 
of  seals  are  hauled-out  on  land,  In  the  spring,  the  number  of  seals  on  land  increases, 
and  wc  can  observe  outbreaks  of  various  final  sizes.  In  the  summer,  the  majority  of  the 
population  is  hauled  out,  the  final  size  is  determined  by  the  law  of  large  numbers  and 
wre  can  observe  major  outbreaks.  The  mean  final  size  of  outbreaks  gradually  decreases 
in  the  fall.  In  short,  for  sandy  haul-out  pattern,  and  parameters  used  in  simulations, 
there  can  be  no  outbreaks  in  the  winter,  whereas  at  other  times  of  the  year  there  can 
be  outbreaks  of  various  final  sizes.  In  the  rocky  regions,  seals  are  hauled-out  in  large 
numbers  throughout  the  year,  so  the  possibility  of  a  large  outbreak  exists  year-round. 

Qualitatively  similar  dynamics  are  observed  for  different  probability  of  transmis¬ 
sion  p  the  final  size  of  infection  is  smaller  in  the  winter  than  in  the  summer,  and 
for  the  same  initial  infection  date  the  size  of  the  outbreak  is  larger  in  the  rocky  than 
in  sandy  regions.  For  smaller  value  of  p  (Figure  3-7)  the  final  size  of  the  outbreaks  in 
the  rocky  regions  drops  in  the  winter,  and  the  envelope  defined  by  the  5th  and  95th 
percentile  becomes  wider. 

It  is  intuitive  that  when  haul-out  behavior  varies  with  season,  the  expression  for 
IZo  must  incorporate  the  timing  of  the  outbreak.  Let  l0  be  the  day  when  the  first 
infectious  seal  is  introduced  to  an  otherwise  susceptible  population  of  size  s(j .  The 
probability  that  this  infectious  seals  is  on  land  on  day  t.Q  is  equal  to  hlo .  and  the 
average  number  of  newr  infections  during  day  l0  is  ps0(hto)2.  The  probability  that 
this  infective  is  on  land  on  its  second  day  of  infection  is  /i(o+I.  The  average  number  of 
new  infections  on  the  second  day  of  the  infection  is  p (hto+i)2.  Taking  the  sum  of  all 
the  new  infections  caused  by  one  infectious  individual  throughout  its  entire  infectious 
period  that  lasts  12  days,  gives 


*0+  1 1 


(3.6) 


that  is  different  for  different  days  of  the  year. 
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Figure  3-8:  7£o(fo)  and  the  timing  of  virus  introduction  (dashed  line).  Sq  =  1500. 
Ph  =  0.001,  rn  =  0.0.  Figure  also  shows  the  mortality  (solid  line)  for  sandy  regions 
from  Figure  3-7  and  the  threshold  level  TZq  —  1  (horizontal  line)  for  reference. 

Figure  3-8  shows  7?(,  against  f0  corresponding  to  the  parameters  of  the  sandy 
areas  shown  in  Figure  3-7.  We  would  expect  that  the  region  of  elevated  mortality 
corresponds  with  the  region  of  graph  where  IZo(to)  >  1.  However,  it  turns  out  that 
TZq  predicts  the  final  mortality  poorly  in  this  case,  and  increase  in  mortality  occurs 
occurs  before  TZq  reaches  the  levels  beyond  tin*  threshold. 

The  graph  3-8  show  the  mortality  at  the  end  of  a  300-day  long  outbreak.  For 
example,  take  day  110  as  the  day  the  virus  is  introduced  to  the  population.  On  day 
110,  I  introduce  one  infected  individual  to  an  otherwise  susceptible  population  and. 
for  a  given  set  of  parameters.  I  simulate  10.000  epidemics  each  lasting  300  days.  The 
mortality  observed  at  the  end  of  those  300  days  is  plotted  in  the  graph  with  110  on 
the  x-axis. 

The  value  of  7?o(/o)  is  calculated  based  on  the  haul-out  levels  during  the  infectious 
period;  77o(110)  only  accounts  the  haul-nut  levels  during  days  110-121.  In  the  Fig¬ 
ure  3-5.  we  can  set'  that  if  the  virus  is  introduced  on  day  110,  some  of  the  trajectories 
will  result  in  an  outbreak.  However,  the  number  of  dead  in  these  trajectories  does 
not  begin  to  increase  until  100  days  after  the  beginning  of  the  outbreak  (whereas  the 
number  of  dead  increases  much  sooner  for  days  of  introduction  100  and  200).  This 
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late  increase  in  number  dead  reflects  the  elevated  fraction  of  the  population  hauled- 
out  well  after  the  introduction  day.  If  the  number  of  seals  on  land  is  low.  the  infection 
can  linger  in  the  population  at  a  low  level,  and  turn  into  an  outbreak  well  after  the 
initial  infection  day.  Since  the  value  of  7£0(/u)  only  incorporates  the  changes  in  ht  for 
the  duration  of  one  infectious  period  (12  days),  it  does  not  refleet  the  ‘delay’  in  the 
increase  in  mortality. 

3.3  Estimation  of  1Z{)  from  data 

When  the  initial  population  size  (sQ)  the  haul-out  behavior  and  the  date  of  the  first 
infection  are  known,  estimating  7 reduces  to  estimating  the  probability  of  infection 
p  The  use  of  maximum  likelihood  methods  for  estimating  p  would  require  knowledge 
of  the  number  of  seals  in  each  class  for  every  day  of  the  epizootic.  Let  /()  be  defined 
as  above  (day  of  virus  introduction),  and  let  I  indicate  the  duration  of  the  outbreak 
in  days.  The  likelihood  of  an  epidemic  trajectory  is 


(3.7) 


where  f  is  the  binomial  probability  density  function. 


Had  we  observed  the  number  of  seals  in  each  epidemic  compartment  throughout  the 
outbreak,  and  the  number  of  infectious  and  susceptible  seals  on  land  throughout  the 
outbreak,  we  could  calculate  L(p)  and  estimate  the  value  of  p.  as  the  value  f>  that 
maximizes  this  likelihood. 

The  only  observation  of  the  epidemic  data  in  case  of  PDV  is  the  information  on 
the  stranded  carcasses.  The  number  of  seals  that  gets  stranded  and  reported  will 
depend  on  many  factors  such  as  the  weather  conditions  and  reporting  effort,  so  it 
will  include  an  observation  error.  The  cumulative  number  of  stranded  seals  provides 
the  total  number  of  seals  that  have  died  of  distemper  in  a  particular  location.  A 
more  reliable  way  to  obtain  the  total  mortality  is  from  the  difference  in  census  data 
before  and  after  the  outbreak.  As  in  Chapter  2.  I  let  the  observation  error  be  equal 
to  the  ratio  of  the  total  number  recovered  stranded  seals  to  the  difference  in  census 
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data.  I  assume  the  observation  error  is  constant  throughout  the  outbreak  and  scale 
the  epidemic  curve  to  match  the  total  number  of  seals  that  died  according  to  the 
population  counts.  The  scaled  daily  counts  of  dead  seals  are  then  used  to  construct 
the  estimates  of  xt  and  it. 

To  reconstruct  the  series  of  incidence.  I  equated  the  observed  mortality  with  its 
expectation  under  (2.2(3)  and  find 


xt  = 


Vt±  is 

rn 


(3.9) 


The  incidence  is  integer- valued.  Rounding  x  to  the  nearest  integer  introduces  the 
possibility  that  the  number  of  the  total  individuals  infected  is  larger  than  the  initial 
susceptible  population  size.  Therefore,  to  keep  x  series  in  integer  form.  I  round  the 
right-hand  side  of  the  equation  (2.30)  to  the  nearest  integer  towards  minus  infinity 
using  the  Matlab  command  floor (). 

After  obtaining  xt.  I  use  assumptions  of  the  model  (3,1)  to  reconstruct  the  series 
for  susceptible  seals 

st+\=st  —  xt,  §i  =  N.  (3.10) 

and  the  series  of  infectious  seals  via 


H  —  h  t- 1  — 


!Jt  |  12  -  Vt  , 
m 


ir  =  0. 


(3.11) 


If  the  haul-out  pattern  h,  is  known,  we  can  use  it  together  with  s  and  /  to  estimate 
^  and  I  bv  taking  the  expected  value  of  (3.2), 


.h,  =  f  loor(/r(  .Sf), 
U  =  floor  (ht.it) 


(3.12a) 

(3.12b) 


We  then  treat  estimates  s.x.i  as  though  they  were  actual  observations  and  esti¬ 
mated  p  by  maximizing  t  lie  log-likelihood  of  the  estimates 


I 

f(p)  =  ln  [f  1  -  (i 


k  0 


(3.13) 


over  p. 
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3,3.1  Accuracy  and  the  precision  of  estimation 

Figures  3-9  3-11  summarize  the  accuracy,  precision,  and  the  bias  of  the  estimates, 
for  constant  haul-out  behavior  and  two  different  patterns  of  seasonal  haul-out.  In 
these  figures,  each  box  plot  represents  a  summary  of  7?0  estimates  from  1.000  simu¬ 
lated  epidemic  trajectories  with  known  parameters.  The  value  of  IZq  used  in  simulated 
trajectories  is  indicated  by  a  horizontal  black  line.  When  TZq  >  1.  the  distribution  of 
final  sizes  is  bimodal,  consisting  of  non-outbreaks  and  outbreaks.  To  unambiguously 
distinguish  between  an  outbreak  and  a  non-outbreak,  I  set  a  threshold  of  what  con¬ 
stitutes  an  outbreak  to  be  20%  of  the  expected  deterministic  final  size  1  —  cxp(— Ro). 
and  discard  any  trajectories  that  did  not  reach  this  threshold. 

As  observed  in  the  previous  Chapter,  the  method  is  negatively  biased  When 
haul-out  behavior  is  present  this  bias  is  stronger,  especially  for  small  population  size 
(100  individuals)  and  for  small  haul-out  levels.  For  constant  haul-out,  estimates  are 
more  accurate,  precise,  and  less  biased  for  larger  population  sizes.  When  the  haul- 
nut  behavior  varies  in  time,  the  negative  Idas  remains  strong  for  the  population  sizes 
observed  in  the  data  set. 

The  main  source  of  error  is  again  a  numerical  one.  In  Chapter  2  main  source  of 
bias  was  in  the  round-off  error  in  the  reconstruction  of  the  incidence  data.  In  addition 
to  underestimating  incidence  by  rounding  off  the  equation  (3.9)  to  the  nearest  integer 
towards  zero,  there  are  two  other  round-off  steps  in  equations  (3.12)  that  further 
underestimate  the  transmission. 

3.4  Application  to  phocine  distemper  data 

After  developing  the  model  with  the  haul-out  behavior,  studying  its  dynamics,  and 
deriving  a  method  to  estimate  its  epidemiological  parameters,  I  want  to  apply  this 
methodology  to  the  available  phocine  distemper  virus  data.  Here  1  focus  on  those  seal 
populations  for  which  both  the  haul-out  data  and  the  epidemic  curves  an*  available 
for  both  1988  and  2002  outbreaks:  Anholt,  German  Wadden  Sea  (WSNS),  Dutch 
Wadden  Sea  (WSNL),  and  Moray  Firth.  The  population  sizes  before  the  outbreaks 
fur  these  locations,  the  total  number  of  seals  that  have  died,  and  initial  infection  dates 
are  listed  in  Table  3.1. 

The  information  on  the  probability  of  being  on  land,  ht.  is  obtained  from  the 
studies  of  the  haul-out  behavior  of  seals. 
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Figure  3-9:  Estimation  accuracy  of  the  method  for  estimating  TC(|  from  epidemic 
curves  of  the  model  with  constant  haul-out.  for  various  values  of  Rq  and  So-  Each 
box  plot  represents  a  summary  of  estimates  from  1.000  simulations  using  the  R{)  value 
indicated  by  the  black  horizontal  line.  The  box  represents  the  inter-quart ile  range  of 
the  estimates,  and  the  red  line  is  the  median.  The  whiskers  are  lint's  extending  from 
each  end  of  the  box  to  show  the  extent  of  the  rest  of  the  data;  the  maximum  length 
of  the  whiskers  is  1.5  times  the  inter-quartile  range.  Data  that  fall  beyond  the  ends 
of  whiskers  are  are  shown  with  red  plus  signs. 
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Figurt*  3-10:  Estimation  accuracy  of  the  method  for  estimating  Tla  from  epidemic 
curves  of  the  model  with  haul-out  pattern  for  sandy  regions,  for  various  values  of  TZ0 
and  Sq.  Each  box  plot  represents  a  summary  of  estimates  from  1,000  simulations 
using  the  TZ$  value  indicated  by  the  black  horizontal  line.  Box  plots  are  constructed 
as  described  in  Figure  3-9. 
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Figure  3-11:  Estimation  accuracy  of  tile  method  for  estimating  7 ?()  from  epidemic 
curves  of  the  model  with  haul-out  pattern  for  rocky  regions,  for  various  values  of  7v(i 
and  S[).  Each  box  plot  represents  a  summary  of  estimates  from  1 ,000  simulations 
using  the  7^n  value  indicated  by  the  black  horizontal  line.  Box  plots  are  constructed 
as  described  in  Figure  3-9. 
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Tabic  3.1:  Sizes  of  harbor  seal  populations  before  1988  and  2002  PDV  outbreaks  for 
the  regions  whose  epidemic  and  haul-out.  curves  were  used  in  estimating  Ho  /J{oc) 
is  the  total  number  of  seals  that  have  died  in  each  location,  and  to  designates  the 
Julian  day  of  15  days  before  the  day  when  the  first  dead  seal  was  recovered  in  each 
location.  Data  provided  by  Karin  Harding  and  Tern  Harkonen. 


N 

1988 
D(o o) 

to 

N 

2002 

/l(ac) 

to 

Anholt 

863 

477 

88 

1467 

200 

109 

WSNS 

4602 

2633 

147 

10042 

4689 

183 

WSNL 

1800 

914 

149 

7002 

3251 

152 

Moray  Firth 

1598 

93 

178 

1198 

86 

238 

3.4.1  Haul-out  behavior 

Harbor  seals  have  been  studied  at  their  haul-out  sites  for  decades.  During  1979- 
1986  aerial  surveys  were  conducted  simultaneously  for  Swedish  and  Danish  haul-out 
locations,  which  were  photographed  in  the  peak  haul-out  season.  Seals  were  later 
counted  from  the  photographs  (Heide- Jorgensen  et  ai.  1992;  Harkonen  et  ai.  1999. 
2002).  The  observations  were  reintroduced  after  the  1988  PDV  outbreak. 

Systematic  observations  and  counts  of  seals  have  been  carried  out  on  Anholt  since 
1978  (Heide- Jorgensen  &  Harkonen.  1988).  Sand  dunes  at  the  beaches  where  the 
seals  haul  out  were  used  as  platforms  for  observations.  Seals  were  counted  twice  on 
every  occasion,  and  the  average  number  was  reported.  When  several  observations 
were  reported  for  the  same  day,  I  consider  only  the  maximum  value  for  that  day. 
Each  year,  the  number  of  seals  hauled-out  peaks;  Figure  3-12  shows  the  hauled-out 
data  as  a  proportion  of  that  maximum. 

From  telemetry  studies  of  seals  equipped  with  VHF  transmitters,  we  can  infer 
that  the  proportion  of  the  total  population  hauled-out  at  maximum  is  between  65- 
71%  excluding  pups  (Thompson  et  ai.  1997;  Ries  et  al..  1998;  Harkonen  et  ai.  1999) 
or  between  57-59%  including  pups  and  assuming  that  pups  of  the  year  haul  out  10% 
of  their  time  (hiring  surveys  (Karin  Harding,  personal  communication). 

The  four  locations  I  am  considering  here,  differ  in  the  number  of  observations  of 
haul-out  behavior.  For  Anholt,  there  are  12  years  worth  of  observations,  for  Dutch 
Wadden  Sea  there  is  only  one  year,  whereas  for  Scotland  and  German  Wad  don  Sea 
the  data  is  in  the  form  of  monthly  averages.  In  order  to  have  the  data  for  all  locations 
in  the  same  form,  I  calculate  the  monthly  averages  for  Anholt  and  Dutch  Wadden 
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Figure  3-12:  Observation  data  for  haul-out  behavior  of  seals  on  Anholt  for  the  years 
1978-1986  and  1989-1991. 
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Figure  3-13:  Haul-out  data  for  German  and  Dutch  Wadden  Sea  and  Scotland 
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Figure  3-14:  Monthly  averages  of  haul-out  data  for  Anholt,  German  Wadden  Sea. 
Dutch  Wadden  Sea  and  Scotland. 


95 


(a)  1988  outbreak 


o 

o> 

% 

g 

Vi 

LU 


6 

4 

2 

0 


An  holt 


WSNS  WSNL 

(b)  2002  outbreak 


Moray  Firth 


Figure  3-15:  Estimation  of  TIq  from  epidemic  curves  and  haul-out  data  for  Anholt. 
Nieder-Sadisen  Wadden  Sea  (WSNS).  Dutch  Wadden  Sea  (WSNL).  and  Moray  Firth. 
Solid  lines  represent  TZlt  values  estimated  from  data,  and  box-plots  illustrate  accuracy 
and  precision  of  the  method,  based  on  1000  epidemic  trajectories  simulated  for  each 
location.  Final  size  assumed  to  be  equal  to  1  in  all  plots. 


Sea.  In  order  to  estimate  the  transmission  probability  for  the  model  with  haul-out 

behavior,  we  need  to  know  hi  for  every  day  of  tiie  year  (/  —  l . 305).  Therefore.  I 

use  linear  interpolation  to  missing  data  points,  and  scale  them  so  that  tin1  obtained 
haul-out  curves  iti  Figure  3-14  peak  according  to  the  telemetry  levels 


3,4.2  T^o  estimates  for  1988  and  2002  PDV  outbreaks 

Using  epidemic  curves  from  Harkbneii  fit  al.  (2000),  population  levels  and  tu  values 
listed  in  Table  3.1,  and  the  haul-out  curves  in  Figure  3-14.  I  estimated  values 
for  1988  and  2002  PDV  outbreaks  for  Anholt,  German  (Nieder-Sachsen  region)  and 
Dutch  Wadden  Sea,  and  Moray  Firth.  The  probability  of  death  m  was  calculated 
assuming  that  all  seals  become  infected  over  the  course  of  the  outbreak  (final  size 
1).  Table  3.2  summarizes  7£u  estimates  with  and  without  accounting  for  the  haul-out 
behavior. 

I  evaluated  the  accuracy,  precision  and  bias  of  each  7Ui  estimate  by  simulating 
1,000  epidemic  trajectories  using  the  values  of  parameters  estimated  for  each  location 
Since  tin1  resulting  trajectories  consist  of  both  non-outbreaks  and  outbreaks.  I  discard 
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Table  3.2:  Comparison  of  the  TZq  estimates  with  and  without  accounting  for  haul-out 
behavior. 


Ill 

1988 

no  haul-out 

haul-out 

m 

2002 

no  haul-out 

haul-out 

Anhoit 

0.55 

2.92 

1.05 

0.14 

2.22 

1.00 

WSNS 

0.57 

2.51 

1.17 

0.47 

2.71 

5.36 

WSNL 

0.52 

2.33 

1.94 

0.41 

2.5 

1.77 

Moray  Firth 

0.06 

2.19 

2.91 

0.07 

1.02 

2.00 

any  trajectories  that  did  not  reach  t lie  20%  of  the  expected  deterministic  final  size 
1  —  exp(— /?o).  Estimates  of  TZu  from  the  remaining  simulated  trajectories  are  shown 
with  box  plots  in  Figure  3-15  for  both  outbreaks. 

The  simulations  suggest  that  bias  is  negative  in  Anholt  and  in  Moray  Firth  for 
both  1988  and  2002.  Further,  the  simulations  suggest  that  TZo  estimates  are  both  more 
accurate  and  more  precise  for  the  Wadden  Sea  populations.  This  may  be  because  the 
population  size  in  the  Wadden  Sea  is  larger  than  in  Anholt  or  Moray  Firth,  especially 
in  2002  when  the  German  Wadden  Sea  population  alone  numbers  over  10,000.  Figures 
3-9- 3- 11  all  indicate  the  estimates  are  more  accurate  and  more  precise  for  large  Sq. 

The  fractions  of  the  population  that  died  in  the  epidemic  (I  refer  to  this  fraction 
as  the  final  morality)  vary  among  locations  in  both  PDV  outbreaks.  Since  each  year 
1’DV  appeared  only  once  in  each  location,  we  do  not  know  whether  the  mortality 
within  tiie  location  would  also  vary,  had  there  been  multiple  outbreaks  throughout 
the  same  year.  To  study  how  mortality  would  vary  within  a  location  for  different  days 
of  virus  introduction,  I  simulated  multiple  epidemic  trajectories  using  the  parameters 
for  Anholt,  WSNS,  WSNL  and  Moray  Firth  corresponding  to  1988  (Figure  3-10) 
and  2002  outbreak  (Figure  3-17).  For  both  outbreaks,  the  observations  of  mortality 
provide  the  probability  of  death  m,  so  the  data  falls  on  the  upper  bound  of  the 
simulated  trajectories. 

In  all  locations,  except  for  Moray  Firth,  in  both  the  1988  and  the  2002  PDV- 
oi  it  break  scenario  the  final  mortality  is  zero  at  the  beginning  of  the  year,  starts  to 
increase  between  days  50  and  100,  peaks  around  day  200.  and  drops  back  to  zero  at 
the  end  of  the  year.  The  peak  in  mortality  occurs  before  the  peak  in  the  haul-out 
behavior.  The  mortality  in  the  Moray  Firth,  although  small,  is  constant  throughout 
the  year.  The  final  sizes  of  the  simulated  outbreaks  start,  also  at  zero  at  the  beginning 
of  the  year,  begin  to  increase  after  day  50.  peak  before  200,  and  then  drop  down  to 
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Figure1  3-16:  Simulations  of  PDV  outbreaks  for  Anholt,  WSNS,  WSNL,  and  Moray 
Firth  using  the  parameters  corresponding  to  the  1988  outbreak  (Tables  3.1  and  3.2). 
For  each  day  of  the  year  I  simulated  1,000  epidemic  trajectories  using  that  day  as  the 
day  of  virus  introduction,  and  summarized  the  final  mortalities  and  final  sizes  with 
the  mean  (solid  line)  and  with  the  5th  and  95th  percentile  (gray  envelope).  Vertical 
line  and  the  circle  represent  15  days  before  the  first  dead  sen!  was  found  and  observed 
final  mortality  in  each  location. 
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Figure  3-17:  Simulations  of  PDV  outbreaks  for  Anholt,  WSNS,  WSNL,  and  Moray 
Firth  using  the  parameters  corresponding  to  the  2002  outbreak  (Tables  3.1  and  3.2). 
For  each  day  of  the  year  1  simulated  1,000  epidemic  trajectories  using  that  day  as  the 
day  of  virus  introduction,  and  summarized  the  final  mortalities  and  final  sizes  with 
the  mean  (solid  line)  and  with  the  5th  and  95th  percentile  (gray  envelope).  Vertical 
line  and  the  circle  represent  15  days  before  the  first  dead  seal  was  found  and  observed 
final  mortality  in  each  location. 
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zero.  But.  in  the  Moray  Firth,  almost  all  values  of  final  sizes  are  possible  throughout 
the  year,  which  is  why  mortality  is  positive  throughout  the  year  in  that,  location. 
Since  the  haul-out  levels  are  well  beyond  zero  throughout  the  year  in  tin;  Moray 
Firth,  there  is  enough  interaction  between  the  seals  to  promote  the  transmission  of 
the  virus  year  round.  The  level  of  mort  ality  in  Moray  Firth  is  low  because  of  the  tow 
probability  of  death,  m,  estimated  for  that,  location.  It  is  only  6%  whereas  m  is  over 
50%  in  other  three  locations  for  the  1988  outbreak. 

In  addition  to  variation  in  the  final  mortality  throughout,  the  year,  Figures  3- 
10  and  3-17  point  out  there  can  be  significant  variation  in  mortality  even  for  the 
same  day  of  virus  introduction.  During  the  first  half  of  t  he  year,  when  the  fraction 
of  the  population  that  is  hauled-out  increases,  the  envelope  determined  by  the  5th 
and  95th  percentiles  is  very  wide  and  includes  outbreaks  of  all  possible  final  sizes. 
When  the  numbers  on  land  are  decreasing,  this  envelope  is  narrow  and  filial  sizes  are 
distributed  around  some  deterministic  value.  Even  though  the  haul-out  curve  has 
a  roughly  symmetric  shape,  the  envelopes  arc  not  symmetric.  This  is  because  the 
total  number  of  infections  does  not  depend  only  on  the  fraction  of  seals  hauled-out 
on  the  initial  infection  day.  but  also  on  the  fraction  of  seals  hauled-out  throughout 
the  outbreak. 


3.5  Discussion 

The  combination  of  the  seasonal  behavior  of  seals  and  the  timing  of  the  virus  in¬ 
troduction  alone  can  explain  the  large  differences  in  mortality  among  regions.  If  the 
virus  is  introduced  to  the  population  in  the  winter  when  the  population  levels  on  land 
are  low,  there  will  be  a  small  outbreak  and  the  population  will  suffer  low  mortality. 
A  large  outbreak  is  most  probable  in  the  summer,  before  the  number  of  the  seals  on 
land  peaks* 

The  importance  if  the  timing  of  the  virus  introduction  and  its  influence  on  tile 
mortality  cannot  be  detected  unless  the  seasonal  behavior  is  present  in  the  model.  The 
importance  of  seasonality  has  been  well  documented  for  other  diseases.  For  childhood 
diseases  like  measles  and  chicken  pox,  seasonality  comes  from  the  aggregation  and 
dispersal  of  schoolchildren  during  and  after  the  school  year  {/..</*,  Anderson,  1996; 
Bjornstad  et  ai ,  2002).  For  vector-borne  diseases  such  as  malaria,  the  seasonality 
comes  from  the  fluctuations  in  the  mosquito,  i *  e,,  vector,  abundance  (Anderson, 
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19%).  In  all  of  these  examples,  the  seasonality  influences  the  transmission  of  the 
disease  and  ignoring  it  can  lead  to  wrong  conclusions  about  the  dynamics  of  the 
epidemic  in  the  population. 

Haul-out  behavior  and  the  timing  of  the  virus  introduction  not  only  explain  the 
variation  in  mortalities  among  regions,  but  can  also  explain  the  differences  in  the 
mortality  between  two  outbreaks  at  the  same  location.  Seasonality  is  not  the  only 
explanation  for  the  difference  in  mortality  among  locations.  Two  other  interpretations 
are:  (i)  Differences  in  mortalities  are  linked  to  pollution,  because  mortality  rates  are 
higher  in  regions  with  higher  concentration  of  PCBs,  pollutants  known  to  suppress 
the  immune  system  of  many  animals  (Bergman  et  al, ,  1992;  Mortensen  et  al. .  1992; 
De  Swart.  1995:  de  Koeijer  ef  al.,  1998).  (ii }  Harbor  seal  populations  are  genetically 
differentiated,  and  different  gene  frequencies  could  lead  to  different  susceptibility  of 
different  sub-populations  and  influence  the  mortality  of  local  populations  (Stanley 
et  al. ,  1996;  Goodman,  1998). 

Even  though  the  combination  of  the  seasonality  of  transmission  and  the  timing 
of  the  virus  introduction  clearly  play  a  substantial  role  in  determining  the  final  size 
and  the  final  mortality  of  an  outbreak,  other  factors  cannot  be  entirely  ruled  out. 
Differences  in  pollution  levels  exist,  anti  many  pollutants  have  proven  innmmotoxic 
effects.  Immunosuppressants,  such  as  PCBs,  and  different  inherent  susceptibility  to 
disease  can  elevate,  or,  in  the  case  of  decreased  susceptibility,  lower  t  lie  levels  of 
mortality  predicted  by  the  model.  However,  I  think  mortality  levels  “correcting” 
for  immunosuppression  and  genetic  differentiation  would  fall  within  or  close  to  the 
bounds  described  by  the  model  with  seasonal  behavior  alone. 

Contamination  with  organochlorines  may.  however,  play  an  important  role  in  de¬ 
termining  the  time  a  certain  population  takes  to  recover  from  such  a  serious  mortality 
event,  since  organochlorine  pollution  can  lower  the  reproductive  success  of  seals  (Rei- 
jnders,  1986,  2003).  Growth  rate  of  harbor  seals  is  already  constrained  by  a  single 
birth  per  female  per  year  (Harkonen  &  Heide- Jorgensen,  1990),  so  any  further  de¬ 
crease  in  growth  rate  due  to  lowered  reproduction  success  can  lead  to  much  slower 
recovery  of  the  population  which  can  be  hazardous  in  the  case  of  recurrent  virus 
outbreaks. 
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Abstract 


The  dynamics  of  HIV/AIDS  epidemics  in  a  certain  region  is  determined  not  only 
by  virology  and  virus  transmission  mechanisms,  but  also  by  region's  socioeconomic 
aspects.  In  this  paper  we  study  the  HIV  transmission  dynamics  for  Cuba.  Wo  mod¬ 
ify  the  model  of  Arazoza  k  Lounes  (2002)  according  to  the  background  about  the 
virology,  as  well  as  the  socioeconomic  factors  that  impact  the  epidemiology  of  the 
Cuban  HIV  outbreak.  The  two  main  methods  for  detection  of  HIV/AIDS  cases  in 
Cuba  are  ‘random’  testing  and  contact  tracing.  As  the  detection  equipment  is  costly 
and  depends  on  biotechnological  advances,  the  testing  rate  can  be  changed  by  many 
external  factors.  Therefore,  our  model  includes  time-dependent,  testing  rates.  By 
comparing  our  model  to  the  1986-2000  Cuban  HIV/AIDS  data  and  de  Arazoza  and 
Lounes  model,  we  show  that  socioeconomic  aspects  are  an  important  factor  in  deter¬ 
mining  the  dynamics  of  the  epidemic. 


4.1  Introduction 

Human  Immunodeficiency  Virus  (HIV)  is  a  global  problem  with  an  estimated  40  mil¬ 
lion  infected  worldwide  (UNAIDS,  2004).  Population  infect  ivity  estimates  range  as 
high  as  8.5V  for  Sub-Saharan  Africa,  and  as  low  as  less  than  0.1  V  for  East  Asia 
and  Australia/New  Zealand.  Cuba,  in  this  respect,  is  remarkable  as  its  infect  iv- 
ity  is  estimated  as  less  than  0.1%  despite  its  status  as  a  relatively  resource-poor 
nation  (Kirkpatrick,  1997:  A  AW,  2005).  The  understanding  of  Cuban  HIV/ AIDS 
infectivitv  dynamics  may  assist  tin*  design  of  preventive  and  reactive  measures  to 
HIV  in  countries  with  high  HIV  prevalence.  This  hypothesis  is  supported  by  Cuba's 
well-developed  health  care  system  despite  its  resource  limitations  (Kirkpatrick,  1997; 
A  AW,  2005). 

The  purpose  of  this  paper  is  to  develop  a  new  model  that  explains  the  dynamics 
of  HIV/AIDS  epidemic  in  Cuba,  focusing  on  the  period  of  1986-2000.  We  built  our 
model  upon  the  work  of  Arazoza  Lounes  (2002)  and  we  confront  both  models  wit  h 
t  he  available  data  (Arazoza  k  Lounes,  2002).  We  begin  wit  h  a  review  of  the  virology 
of  HIV/AIDS  within  the  socioeconomic  framework  of  Cuba  1986-2000  in  Sect  ion  4.2. 
The  formulation  and  brief  analysis  of  the  mathematical  model  follows  in  Sect  ion  4.3, 
as  well  as  the  comparison  of  the  model  with  data.  We  finish  with  a  discussion  in 
Section  4.4. 


106 


Table  4.1:  New  cases  of  HIV,  AIDS,  AIDS-related  deaths  in  Cuba  1986-2000  (Arazoza 
k  Lounes,  2002). 


Year 

IllV-cases 

AIDS-cases 

Death  due  to  AIDS 

1986 

99 

5 

2 

1987 

75 

11 

4 

1988 

93 

14 

6 

1989 

121 

13 

5 

1990 

140 

28 

23 

1991 

183 

37 

17 

1992 

175 

71 

32 

1993 

102 

82 

59 

1994 

122 

1Q2 

62 

1995 

124 

116 

80 

1996 

234 

99 

92 

1997 

363 

129 

99 

1998 

362 

150 

98 

1999 

493 

176 

122 

2000 

545 

251 

142 

4.2  Background 


With  a  total  population  of  11  million,  and  less  than  1000  infected,  Cuba’s  HIV/AIDS 
epidemic  is  a  small  one.  As  part  of  the  HIV/ AIDS  prevention  program,  Cuba  has 
an  active  search  of  seropositives  through  the  sexual  contacts  of  known  HIV-infected 
persons;  t  his  system  is  called  contact  tracing.  Infected  persons  are  also  found  t  hrough 
a  ‘blind’  search  of  blood  donors,  pregnant  women,  persons  with  other  sexually  trans¬ 
mitted  diseases,  etc.  (Arazoza  k  Lounes,  2002).  Both  methods  arc  very  successful 
in  locating  HIV-positive  persons  (Arazoza  k  Lounes,  2002).  The  numbers  of  newly 
diagnosed  HIV  cases,  AIDS  cases,  and  AIDS-related  deaths  per  year  in  Cuba  are  de¬ 
tailed  in  Table  4.2  and  plotted  in  Fig.  4-1  below  Arazoza  k  Lounes  (2002).  However, 
fluctuations  in  these  numbers  are  due  to  both  the  character  of  the  HIV  virus  and  the 
manner  in  which  the  Cuban  population  has  been  monitored  for  its  presence.  A  model 
which  does  not.  distinguish  between  virology,  the  socioeconomic  framework  which  this 
virology  exists  (i.e.,  the  epidemiology),  and  how  this  framework  has  been  observed, 
may  generalize  very  poorly. 
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Figure  4-1:  New  cases  of  HIV.  AIDS,  AIDS-related  deaths  in  Cuba  1986-2000  (Ar 
zoyai  &  Lounes,  2002). 


4.2.1  Virology 


An  average  HIV-infected  individual  progresses  through  distinct  stages  of  the  disease. 
The  infeetiousness  (i.e.,  the  probability  of  transmission)  varies  greatly  depending 
upon  the  stage  of  the  disease.  First  comes  a  period  of  primary  infection  (lasting  part 
of  a  year  Ahlgren  et  ai,  1990).  During  the  primary  stage,  infeetiousness  first  rises 
and  then  drops.  Seroconversion  usually  occurs  before  the  end  of  the  first  year.  HIV  is 
an  asymptomatic  period  (Ahlgren  et  ai,  1990,  averaging  7  years  without  treatment) 
in  which  infeetiousness  is  low.  This  is  followed  by  a  symptomatic  stage  (averaging 
three  years  until  death  without  treatment  Ahlgren  et  ai,  1990)  where  infeetiousness 
rises  again.  Although  toward  the  end  of  the  symptomatic  stage  individuals  are  expe¬ 
riencing  severe  AIDS  and  activity  is  decreased,  the  symptomatic  stage  begins  while 
individuals  are  relatively  healthy  and  still  very  active.  The  average  stage  infectivity 
rat  es  for  semen  has  a  of  a  small  peak  shortly  after  init  ial  infection  followed  by  a  larger 
peak  during  the  symptomatic  phase  (Rapatski  et  ai,  2005).  This  correlates  with  the 
changes  in  viral  load  observed  as  a  person  progresses  through  the  disease  (Pantaleo 
et  ( li,  1993;  Clark  et  ai,  1991;  Darr  et  ai ,  1991;  Anderson,  1990).  This  pattern  is  due 
to  the  physiology  of  the  disease,  the  way  the  infected  persons’  bodies  interact  with  the 
virus  (Gray  et  ai,  2001;  Saracco  et  ai,  1993;  Piatak  et  ai,  1993;  Vincenzi,  1994),  and 
is  largely  independent  of  the  sexual  practices.  In  Cuba,  most  of  the  transmissions 
occur  through  sexual  intercourse  (about  a  1:1  ratio  of  heterosexual  to  homosexual 
transmission,  Holtz,  n.d.;  Hsieh  et  ai,  2001). 


4.2.2  HIV  in  Cuba 

Cuba  treated  the  introduction  of  HIV  into  the  country  in  1986  as  a  public  health 
emergency,  introducing  control  measures  to  contain  the  spread  of  the  disease.  As  a 
result.  Cuba  has  one  of  the  lowest  prevalence  rates  of  HIV  infection  in  the  world. 
Cuba’s  HIV  prevalence  of  0.03%  is  nearly  11  times  lower  than  that  of  the  United 
States  (Perez-Stable,  1991;  Burr,  1997).  In  1986,  Cuba  introduced  a  national  screen¬ 
ing  program.  Cuba  had  a  well-developed  health  care  system  that  assigned  a  primary 
care  physician  to  all  citizens  and  conducted  routine  surveillance  for  infectious  disease 
(Waitzkin  et.  ai,  1997;  Feinsilver,  1989).  To  reduce  the  risk  of  transmission  Cuba 
instituted  numerous  measures,  including  contact,  tracing,  isolation  (quarantine)  of 
HIV-infected  individuals  and  a  total  ban  on  the  import  of  blood  and  blood  byprod- 
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ucts  (Holtz,  n.d.;  Hsieh  et  al..  2004).  Initially,  quarantine  individuals  lived  in  isolation 
in  sanitariums.  By  1993,  patients  could  choose  between  living  within  a  sanitarium 
or  living  at  home.  In  the  sanitariums,  people  are  provided  with  good  meals,  a  par¬ 
tial  salary,  free  medications  and  care  from  physicians  (Santana  et  al.,  1991).  Most 
individuals  could  not  provide  the  care  necessary  for  them  and  therefore  most  choose 
to  live  in  the  sanitariums  (Holtz,  n.d.}.  Once  a  person  is  quarantined,  they  are  no 
longer  a  factor  in  t  he  transmission  of  the  disease.  Contact  tracing  in  Cuba  involves 
the  search  of  HIV-positive  persons  through  the  sexual  contacts  of  known  HIV-infected 
individuals.  This  practice  has  proven  to  be  quite  effective  in  Cuba  (Hsieh  et  al.,  2004). 
Since  a  significant  fraction  of  those  found  to  be  HIV  positive  occur  through  contact 
tracing,  a  model  of  HIV  in  Cuba  must  allow  for  contact  tracing. 


HIV  Data 


To  model  the  Cuban  HIV  epidemic,  one  has  to  acknowledge  contact  tracing  and 
quarantines  as  well  as  any  inconsistencies  with  the  available  data  (Table  4.2.  Figure  4- 
1).  The  first  column  in  Table  4.2  represents  those  individuals  that  tested  positive  for 
HIV  during  that  year;  t  hey  may  have  acquired  the  disease  some  time  before.  The 
number  of  total  II IV'  eases  in  column  one  includes  both  newly  tested  HIV-positives  and 
the  people  in  the  AIDS  stage.  Because  of  this  combination  along  with  the  aggressive 
testing  of  Cuba,  we  believe  the  AIDS  data  (column  3)  to  be  more  reliable  than  the 
HIV  data.  From  Table  4.2,  it  appears  as  though  from  1990-1992  there  was  an  increase 
in  the  number  of  newly  HIV  infected  persons.  This  increase  was  due  to  the  discovery 
and  contact  tracing,  from  approximately  1990  to  1992,  of  a  highly  sexually  active 
group  (de  Arazoza  et  al..  2003).  Because  of  a  United  States  embargo  in  1992,  new 
HIV  testing  equipment  was  no  longer  available  to  Cuba  (Holtz,  n.d.),  leading  to  a 
decrease  in  the  number  of  newly  HIV  infected  individuals  being  discovered  I  hat.  year. 
These  two  events  are  highlighted  in  Fig.  4-2  A  model  of  the  HIV  epidemic  in  Cuba 
must  account  for  these  two  significant  events. 
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Figure  4-2:  Socioeconomic  factors  that  influenced  diagnosis  of  HIV  positive  persons. 

4.3  Mathematical  models  and  analysis 

4.3.1  Previous  Model 

De  Arazoza  and  Lounes  have  modeled  Cuba’s  I II V/ AIDS  epidemic.  They  consider 
three  divisions  of  the  population,  undiagnosed  HIV  positive  {{/),  diagnosed  HIV  posi¬ 
tive'  ( D }.  and  AIDS  (A)  with  the  following  constant  coefficients  (values  listed  in  Table 


1.  N,  total  size  of  the  sexually-active  population, 

2.  at.  the  rate  of  recruitment  of  new  HIV- infected  persons,  infected  by  U . 

3.  o',  t  he  rate  of  recruitment  of  new  HIV-infected  persons,  infected  by  D. 

4.  k i,  the  rate  at  which  the  unknown  HIV-infected  persons  are  detected  by  the 
system  (“random”  search), 

5.  fc2,  the  rate  at  which  the  unknown  HIV-infected  persons  are  detected  through 
contact  tracing, 
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6.  /?,  the  rate  at  which  the  HIV  positives  develop  AIDS, 

7.  ft,  the  mortality  rate  of  the  sexually  active  population, 

8.  ft',  the  mortality  rate  of  the  population  with  AIDS. 

In  this  model  there  are  two  ways  individuals  can  go  from  unknown  HIV  infected 
([/)  to  diagnosed  HIV-infected  ( D).  through  contact  tracing  (k2UD) and  detection 
through  ail  other  random  searching  for  seropositives  {k\U).  Authors  assume  that  the 
known  HIV  infected  persons  arc  infectious,  but  at  a  much  lower  rate  than  those  that 
do  not  know  they  are  infected. 

Their  model  equations  are: 


U'  =  aNU  +  a'ND  -  (k ,  +  ft  +  ff)U  -  hVD. 
D'  =  k)U  +  k2UD-{fi  +  /3)D, 

A '  =  !3{U  +  D)  -  ft' A. 


(4.1a) 

(4.1b) 

(4.1c) 


4,3,2  Model  Design 

To  improve  upon  the  previous  model  by  de  Arazoza  and  Louues,  we  have  made  three 
major  changes: 

1.  We  consider  four  divisions  of  the  population,  susceptible  (5),  undiagnosed  HIV 
positive  {£/),  diagnosed  HIV  positive  { D ),  and  AIDS  (.4).  We  considered  this 
to  be  a  closed  population  and  all  births  equal  deaths. 

2.  We  incorporate  the  variation  in  infect. ivity  as  a  person  progresses  through  the 
disease,  by  considering  the  rate  for  a  susceptible  to  be  infected  by  an  individual 
with  AIDS.  lj.  With  the  aggressive  “random'’  testing  in  Cuba,  by  tin*  time 
individuals  progress  to  the  AIDS  stage  they  have  been  diagnosed.  Although 
individuals  with  AIDS  are  much  more  infectious  than  individuals  with  HIV 
(Rapatski  et  al. .  2005),  an  AIDS  individual  would  have  fewer  contacts  wit  h 
susceptible  persons  compared  to  the  contacts  made  by  undiagnosed  individuals 
with  susceptibles  thus,  u,’  is  lower  than  the  rate  for  a  susceptible  to  be  infected 
by  an  undiagnosed  HIV  person,  denoted  o'.  When  comparing  persons  in  the  the 
AIDS  stage  and  diagnosed  HIV  persons,  since  AIDS  stage  is  more  infectious, 
we  assume  t  he  rate  for  a  susceptible  to  be  infected  by  an  individual  with  AIDS, 
u.'.  to  be  higher  than  the  rate  of  a  diagnosed  HIV  posit  ive  individual,  o'. 
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3.  Undiagnosed  individuals  are  diagnosed  by  their  doctors  at  a  rate  k\ ,  and  through 
cont  act  tracing  at  a  rate  k2.  We  consider  three  phases  for  contact  tracing,  1986- 
1989,  1990-1991  and  1992-2000,  and  two  phases  for  “random”  testing,  1986-1991 
and  1992-2000.  In  each  period,  k\  and  k2  are  constant.  We  estimate  that  during 
1990  and  1991  contact  tracing  increased  25%  because  of  the  detection  of  a  highly 
sexually  active  group,  and  that  after  1992  diagnosis  by  doctors  was  reduced  to 
75%'  of  its  former  value  due  to  the  United  States  embargo.  We  obtain  values  of 
k\  and  k2  by  fitting  the  Cuban  HIV/ A  IDS  data. 

The  dynamics  of  the  Cuban  HIV/AIDS  epidemic  are  described  by  the  following 
model: 


S'  = 

~{uA  +  aU  +  a'D)S  +  fi  A  +  tt{U  +  D), 

(4.2a) 

V*  = 

(uM  +  otV  +  a'D)S 

-  {kl{t)  +  (i  +  0)U  -k2{t)UD, 

(4.2b) 

D'  = 

k\{t)U  +  k2(t)UD  - 

{n  +  P)D, 

(4.2c) 

A'  = 

0{U  +  D)  -  ^ A. 

(4. 2d) 

This  model  holds  within  each  of  the  periods.  The  initial  conditions  for  each  period 
are  taken  to  keep  the  overall  solution  continuous  (i.e.,  initial  conditions  art*  the  final 
conditions  for  the  previous  period).  Solutions  to  (4.2)  with  positive  initial  condi¬ 
tions  remain  positive  for  all  periods.  System  (4.2)  has  a  unique  solution  with  initial 
conditions  (S(U).  t/(0),  £>(0),  A(0))  =  (5.5  million,  230,  94,  3). 


4.3.3  Numerical  Results 


Estimates  of  k\  (t)  and  k2{t)  are  obtained  by  minimizing  the  following  error  function. 
For  each  of  the  fifteen  years  we  compute  the  square  of  the  difference  between  our 
model  epidemic  and  the  Cuban  HIV  data  given  in  Table  4.2.  Let  RMS  denote  the 
square  root  of  the  average  of  those  fifteen  numbers, 


RMS  Error 


I  r  ^  [-f^r nodel(t)  ^UlVDataiJ )] 

.  °  1986-2000 


(4.3) 


We  select  the  values  of  ki(t)  and  k2(t)  that  minimize  RMS,  by  taking  the  gradient 
of  the  {RMS  Error)2  and  using  Newton’s  method  to  find  a  zero  of  the  vector  field. 
The  parameter  values  are  given  in  Table  4.2.  The  initial  values  for  U,  D  and  A 
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Table  4.2:  Values  of  parameters  used  in  simulations. 


Parameter 

Description 

de  Arazoza 

Ours  (Source) 

5(0) 

Initial  condition  for  Sus- 
ceptibles 

N/A 

5.5  million  (a) 

1/(0) 

Initial  condition  for  HIV 
Undiagnosed 

230 

230  (b) 

0(0) 

Initial  condition  for  HIV 
Diagnosed 

94 

94  (b) 

.4(0) 

Initial  condition  for  AIDS 

3 

3  (b) 

Rate  for  a  susceptible  in¬ 
dividual  to  become  in¬ 
fected  by  an  individual 
with  AIDS 

N/A 

8.5  •  10"8  (c) 

a 

Rate  for  a  susceptible  indi¬ 
vidual  to  become  infected 
by  an  undiagnosed  HlV-f 
individual 

9.3267  *  10"8 

9.3267  ■  R)  8  (b) 

o' 

Rate  for  a  susceptible  indi¬ 
vidual  to  become  infected 
by  an  diagnosed  HIV+  in¬ 

5.4  •  nr‘J 

5.4  ■  10”9  (b) 

0 

dividual 

Rate  at  which  HIV+  indi¬ 
viduals  develop  AIDS 

0.10788 

0.14  (d) 

{• 

Mortality  rate  for  HIV  pos¬ 
itive  individuals 

0.75 

0.75  (b) 

l>' 

Mortality  rate  for  individ¬ 
uals  with  AIDS 

0.0053 

0.0053  (b) 

A :i«) 

‘Random’  testing  rate  per¬ 

0.3743 

1986-1991  0.3850  (e) 

formed  by  doctors 

1992-2000  0.2929  (e) 

M0 

Testing  rate  due  to  contact 

2.27-10-5 

1986-1989,  3.26- 10“5 

(e> 

tracing 

1992-2000 

1990-1991  5.89-1 0~*‘ 

(e) 

(a)  UN  (2005) 

(b)  Arazoza  &:  Lounes  (2002) 

(c)  Estimate  based  on  a  and  o' 

(d)  Ahlgren  et  al.  (1090) 

(e)  Estimate  to  fit  data 


Figure  4-3;  Comparison  of  model  (4.2),  and  Arazoza  &  Lounes  (2002)  model  with 
data  for  HIV  positive  cases  in  Cuba. 

were  chosen  to  be  the  same  as  those  used  in  de  Arazoza  and  Lounes  {Arazoza  & 
Lounes,  2002)  and  5(0)  was  estimated  to  be  5.5  million  (assuming  half  of  the  11 
million  population  (UN,  2005)  are  of  a  sexually  active  age).  The  resulting  curves  for 
both  the  diagnosed  HIV  cases  and  AIDS  cases  is  shown  in  Fig.  4-3.  We  compared 
our  model  results  with  de  Arazoza  and  Lounes  model.  As  seen  in  Fig.  4-3,  our  model 
is  a  better  fit.  to  the  data. 


4.3.4  Basic  reproduction  ratio 

The  basic  reproduction  ratio,  7?.(l.  is  a  dimensionless  parameter  that  gives  t  he  expected 
number  of  secondary  eases  per  primary  ease  of  infection  in  an  entirely  susceptible 
population.  As  a  result,  TZo  has  a  threshold  value  equal  to  one,  i.e.,  infection  will 
spread  and  result  in  epidemic  if  TZq  >  1,  whereas  the  infection  will  die  out  if  TZq  <  1. 

Model  (4.2)  has  a  disease  free  equilibrium  (DFE),  £o,  given  by 

ro  :  (5,  U,  D,  A)  =  (50, 0, 0, 0).  (4.4) 

7v0  is  calculated  for  constant  values  of  hi  and  k->,  that  is,  there  is  an  7vn  for  each 
time  period.  We  are  interested  in  looking  at  the  stability  of  a  simpler  model  of  (4.2) 
with  each  k  constant  throughout.  We  are  interested  in  the  final  period,  with  each  k 
set  to  their  final  value. 

From  Diekmann  et  al.  (1990).  lZa  is  the  spectral  radius  (p)  of  the  next  generation 
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matrix  (set1  also  van  den  Driessche  k  Watmongh,  2002),  K. 


TZn  =  p(K), 


(4.5) 


where  K  =  FV"1.  F  and  V  come  from  the  Jacobian  matrix  of  the  linearization  of 
(4.2)  about  tilt1  DFE.  Here,  non-uegai ive  matrix  F  shows  new  infections,  and  the 
inverse  of  the  non- singular  matrix  V  gives  the  expected  t  imes  that  individuals  spend 
in  each  of  the  compartments.  F  and  V  are  respectively  given  by 


ab  a 

fS  u)S ^ 

k\  +  (1  0 

F  - 

0  0 

and  V  = 

-kt-k2D 

1 1+3  0 

U 

0  t)  J 

\  -o 

-a  n'J 

The  basic  reproduction  ration  for  model  (4.2)  is  then  given  by 

'V  _  g(°)  ( n (1  +  /')  +  0w\ 

V°  0  +  p\  k\  +  &+0  it’)' 


(4.G) 


(4-7) 


An  advantage  of  considering  R.t)  on  a  generation  basis,  is  that  we  obt  ain  expression 
(4.7)  for  7vi,  in  terms  of  parameters  of  the  model,  which  provides  implications  for  the 
control  of  the  epidemic  which  we  discuss  in  Section  4.4. 

Since  we  are  int  erested  in  the  simpler  model  where  fe’s  are  const  ant  throughout .  our 
system  becomes  an  autonomous  system.  The  equilibrium  £q  is  locally  asymptotically 
stable  if  7vo  <  1  (van  den  Driessche  Watmongh.  2002),  and  t  he  population  is  not 
vulnerable  to  the  outbreak  of  the  disease.  In  the  case  when  7vn  >  1,  the  DFE  is 
unstable  so  the  disease  can  invade  the  population,  eventually  leading  to  an  endemic 
equilibrium.  These  two  types  of  dynamics  are  illustrated  by  simulations  of  model 
(4.2)  in  Figure  4-4. 

4.4  Discussion 


In  this  paper  we  present  a  new  model  for  studying  H1V/A1DS  epidemic  in  Cuba, 
based  on  the  de  Arazoza  and  bonnes  model  (4.1).  We  modified  their  model  in  three 
ways.  First.,  we  allow  for  “random’’  test  ing  rate  ( k\ )  and  contact  tracing  rate  ( hi }  to 
vary  in  t  ime,  in  order  to  reflect  the  fluctuating  socioeconomic  situation  in  the  coun¬ 
try.  Second,  we  assume  that  persons  who  developed  AIDS  can  infect  the  susceptible 
individuals.  Even  though  the  people  in  the  AIDS  class  have  fewer  sexual  contacts 
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Figure  4-4:  A)  Ko  >  1,  the  system  reaches  endemic  equilibrium,  B)  TZu  >  1.  disease 
free  equilibrium. 
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than  the  asymptomatic,  HIV-positive  individuals,  symptomatic  individuals  are  highly 
infectious.  The  viral  load  in  the  symptomatic  (AIDS)  stage  can  be  up  to  150  times 
higher  than  in  the  asymptomatic  stage  (Rapatski  et  al. ,  2005),  so  the  probability  of 
transmission  of  HIV  remains  substantial  in  the  symptomatic  stage  and  we  include  it 
in  the  model  (parameter  lj).  Lastly,  since  total  population  in  Cuba  is  much  greater 
(more  than  four  orders  of  magnitude)  than  the  number  of  people  affected  by  HIV  and 
AIDS,  Arazqza  &  Lounes  (2002)  assume  that  the  susceptible  population  is  constant 
in  time,  and  thereby  reduce  a  dimension  in  their  system.  We.  on  the  other  hand, 
model  the  changes  in  t  he  susceptible  class  as  well. 

To  test  our  model  we  have  used  the  yearly  HIV-positive,  AIDS  cases,  and  deaths 
due  to  AIDS  in  Cuba  in  the  period  1986-2000  (Table  I  from  Arazoza  &  Lounes  (2002)). 
Data  includes  newly  HIV-infected  people,  t  he  number  of  people  who  developed  AIDS 
symptoms,  and  the  number  of  people  who  died  from  complicat  ions  of  AIDS.  From 
the  data  wo  cannot  infer  the  time  of  HIV-infection. 

The  current  state  of  the  HIV/AIDS  epidemic  in  Cuba  is  described  with  the  pa¬ 
rameter  values  given  in  Table  4.2.  For  these  values,  Km  >  1,  so  the  number  of  new. 
diagnosed  and  undiagnosed,  HIV  infections  in  Cuba  is  increasing.  However,  compared 
with  Ko  values  for  sub  Saharan  Africa  (9.62  Rapatski  et  al.  submitted),  and  India 
(31  Rapatski  et  al. ,  submitted),  Kq  for  the  Cuban  epidemic  is  very  small. 

As  long  as  11m  remains  greater  than  one.  the  IIIV/AIDS  epidemic  will  continue 
to  spread  in  Cuba.  Mechanisms  that  decrease  the  value  of  Km  in  (4.7)  below  the 
threshold  are  the  mechanisms  that  can  put  the  epidemic  under  control.  Equation 
(4.7)  suggest  two  different  ways  of  controlling  the  epidemic:  increasing  the  rate  at 
which  unknown  HIV-infected  persons  are  detected  (k\ ),  and  decreasing  the  rate  of 
infection  (o). 

Let  us  look  at  the  two  possible  mechanisms  of  control  more  closely.  Increasing 
the  testing  rate  requires  more  effective,  precise  and  affordable  HIV-deteetion  tests, 
and  a  thorough  and  systematic  testing  organized  bv  the  public  health  system.  As 
increasing  the  detection  rate  depends  on  advances  in  biotechnology  and  the  structure 
of  the  public  health  system,  increasing  testing  rate  enough  to  bring  Km  below  the 
threshold  is  unlikely  at  the  moment.  On  the  other  hand,  there  are  widely-available. 
affordable  methods  that  decrease,  the  infection  rate,  a.  The  proper  usage  of  condoms 
has  been  shown  to  reduce  the  risk  of  transmission  of  HIV  in  two  ways.  Condom 
usage  reduces  the  risk  of  transmission  of  HIV  itself,  but  it  also  significantly  reduces 
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the  risk  of  transmission  of  other  sexually  transmitted  infections  (STls),  Since  many 
STIs  can  cause  abrasion  of  the  genital  skin  and  membranes,  STIs  may  facilitate  both 
transmission  and  acquisition  of  HIV  (Moss  et  ai,  1987).  Given  that  less  than  a 
third  of  people  use  condoms  with  their  non-regular  partners  (Gardner  et  al..  1999), 
increased  condom  usage  is  a  promising  measure  against  future  spread  of  HIV/ A  IDS 
epidemic  in  Cuba. 
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Contribution:  Analysis  of  the  Lotfea-Volterra  model  with  discrete  delay. 
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Abstract 


Time  delays  produced  by  dispersal  are  shown  to  stabilize  Lotka- Volterra  predator- 
prey  models.  The  models  are  formulated  as  integrodifferential  equations  that  describe 
local  predator- prey  dynamics  and  either  intrapatch  or  interpatdi  dispersal.  Dispersing 
individuals  may  (or  may  not)  differ  in  the  duration  of  their  trip:  these  differences  art' 
captured  via  a  distributed  delay  in  the  models.  Our  result  s  include  those  of  previous 
studies  as  special  cases,  and  show  that  the  stabilizing  effect  continues  to  operate  when 
the  dispersal  process  is  modelled  more  realistically. 


5.1  Introduction 

Interest  in  the  stability  of  predator-prey  and  host -parasite id  systems  has  continued 
unabated  since  the  theoretical  work  of  Lotka  (192(>),  Volterra  (1931),  and  Nicholson 
and  Bailey  (1935)  and  the  experimental  work  of  Cause  (1934).  The  central  question 
raised  bv  their  work  is  this:  how  do  predator-prey  systems  apparently  persist  stably 
in  nature  when  the  most  basic  models  and  experiments  predict  instability?  The  an¬ 
swer  most  often  given  is  t hat  the  models  and  experiments  omit  processes  that  affect 
stability  in  natural  systems.  To  support  this  answer,  theoreticians  and  experimental¬ 
ists  have  proceeded  to  investigate  the  stability  mediating  effects  of  a  long  list  of  such 
processes  (for  examples  see  May  1973,  Hassell  1973.  Crawley  1992.  and  Mueller  and 
Josh!  2000). 

The  basic  theoretical  tool  in  these  investigations  is  the  system  of  Lotka- Volterra 
equations  for  a  prey  with  populat  ion  density  N(T)  and  a  predator  with  population 
density  P{T): 

^  =  (R-AP)N,  (5,1a) 

”  =  (BN-M)P,  (5.11)) 

In  the  absence  of  predators,  the  prey  population  grows  exponentially  at  the  rate  /?, 
and  in  the  absence  of  prey,  the  predator  population  decays  exponentially  at  the  rate 
j XI.  The  predator-prey  interaction  is  captured  by  linear  functional  and  numerical 

responses,  scaled  by  the  parameters  A  and  B.  The  parameters  B,  A.  B.  and  A/  are 

assumed  to  be  positive. 

The  Lotka- Volterra  predator-prey  model  is  often  criticized  because  its  single,  pos¬ 
itive,  equilibrium  point  is  a  center,  i.e.,  a  “neutrally  stable"  equilibrium  surrounded 
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by  a  family  of  periodic  orbits  whose  amplitudes  depend  on  tiie  initial  population 
sizes.  The  slightest  change  to  the  model’s  structure  typically  results  in  qualitatively 
different  behavior.  For  example,  if  /?  decreases  linearly  with  prey  density  the  equilib¬ 
rium  point  is  stable;  on  the  other  hand,  introducing  a  saturating  (Type  II)  functional 
response  turns  the  equilibrium  into  an  unstable  spiral  point  (Gotelli  1995).  This 
structural  instability,  the  critics  argue,  means  that  the  model  cannot  make  any  pre¬ 
dictions  that  are  robust  enough  to  be  tested.  After  all,  we  know  that  model  (5.1) 
does  not  adequately  describe  even  the  most  highly-controlled  experiments. 

Structural  instability  can,  however,  be  used  to  our  advantage.  In  effect,  it  allows 
us  to  use  the  Lotka-Volterra  model  as  an  exquisitely  sensitive  balance,  wit  h  which  we 
can  determine  the  effects  of  the  processes  that  it  ignores.  So,  when  we  say  that  a  Type 
II  functional  response  is  destabilizing,  we  mean  that  it  destabilizes  the  equilibrium 
point  in  model  (5.1).  Similarly,  when  we  say  that  the  presence  of  carrying  capacity 
for  the  prey  tends  to  be  stabilizing,  we  mean  that  it  stabilizes  the  equilibrium  point. 
There  is  a  long  tradition  of  using  the  Lotka-Volterra  equations  in  this  way  (Murdoch 
and  Oaten  1975),  and  we  continue  that  tradition  here. 

Among  the  many  processes  that  the  Lotka-Volterra  equations  ignore,  those  with 
a  spatial  component  have  always  attracted  attention  (Mueller  and  Joshi  2000).  In 
particular,  the  presence  of  a  metapopulation  structure  (i.e.,  locally  interacting  popu¬ 
lations  coupled  via  dispersal)  can  have  interesting  and  variable  effects  (Hanski  1999). 
Taylor  (1990)  and  Mueller  and  .loshi  (2000)  briefly  review  this  topic. 

The  simplest  metapopulation  model  consists  of  two  habitat  patches.  A  simple 
two-patch  extension  of  the  Lotka-Volterra  model  (5.1)  is  given  by: 

dNi/dT  =  (J?  -  APt)  Ni  +  Dn[Ns  -  Ni],  (5,2a) 

dPi/dT  =  (BNi-  M)  +  Dp[Pj  -  P*],  (5.2b) 

for  i  =1.2  and  j  ^  i  (Comins  and  Blatt  1974).  The  subscripts  indicate  the  patch 
number;  Ds  and  DP  are  the  prey  and  predator  emigration  rates. 

Because  predators  are  often  more  mobile  than  their  prey,  many  authors  have 
studied  a  simplified  version  of  model  (5.2)  with  Dt\  —  0.  Jansen  and  de  Roos  (2000) 
provide  a  concise  review  of  the  dynamics  of  this  model  (also  see  Murdoch  and  Oaten 
1975,  Murdoch  et  al.  1992,  Nisbet  et  al.  1992,  Jansen  1995).  When  DP  >  0,  there  is  a 
constant  per  capita  predator  migration  rate  between  the  patches,  This  coupling  does 
not  change  the  equilibrium  values;  there  is  a  spatially  homogeneous  equilibrium  with 
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population  sizes  in  each  patch  equal  to  their  sizes  in  the  uncoupled  case.  Furthermore, 
coupling  the  two  populations  via  predator  dispersal  does  not  change  the  stability  of 
the  equilibrium  point.  The  equilibrium  is  surrounded  by  a  planar  family  of  unstable 
periodic  solutions  on  the  subspace  defined  by  N ,  —  N?  and  P\  =  P>.  Any  initial 
differences  in  prey  or  predator  population  sizes  between  the  two  patches  eventually 
disappear,  so  orbits  converge  to  this  plane.  Large  amplitude  cycles  in  this  plane  arc 
unstable  to  perturbations  off  the  plane,  while  small  amplitude  orbits  are  stable  to 
off-plane  perturbations.  There  appear  to  be  hcteroelinic  orbits  connecting  the  large 
amplitude  solutions  to  the  small  amplitude  solutions.  Thus  perturbations  to  periodic 
orbits  tend  to  result  in  periodic  orbits  of  smaller  amplitude.  Only  in  this  weak  sense 
('an  predator  dispersal  (as  described  in  model  (5.2)  with  DN  =  0)  be  thought  of  as 
stabilizing.  None  of  the  periodic  orbits  is  asymptotically  stable  (perturbations  in  ihe 
plane  do  not  decay)  and  unless  the  initial  condition  is  set  exactly  at  the  equilibrium 
value  in  each  patch,  the  populations  will  ultimately  cycle  in  synchrony. 

So,  predator  dispersal  by  itself  seems  to  be  insufficient  to  stabilize  the  Lotka- 
Volterra  predat  or- prey  interaction.  But  the  description  of  dispersal  in  model  (5.2)  is 
artificial  in  an  important  way:  dispersers  leaving  one  patch  immediately  appear  in 
the  other  patch.  In  nature,  dispersers  take  a  finite  amount  of  t  ime  to  complete  their 
trip.  During  this  time,  migrating  individuals  are  typically  not  participating  in  the 
predator-prev  interaction  because  the  two  species  are  in  different  places  ( Weiss rr  and 
Hassell  1996,  Weisser  et.  al.  1997). 

In  this  paper,  we  develop  a  general  way  to  explicitly  account  for  individual  travel 
times,  and  show  that  dispersal  is  almost  always  stabilizing  when  an  explicit  travel- 
time  is  incorporated  in  the  model.  We  are  not  the  first  to  demonstrate  this  effect. 
Holt  (1984)  and  Weisser  and  Hassell  (1996)  studied  the  effect  of  dispersal  on  the 
stability  of  a  predator-prey  system  in  a  single  patch.  They  coupled  this  patch  to 
itself  via  constant  per  capita  emigration  {at  rate  E)  and  immigration  {at  rate  I)  into 
and  out  of  a  pool  of  dispersers  {with  density  Q(T)).  When  predators  disperse,  the 
model  has  the  form 

=  (R-AP)N,  (5.3a) 

=  {BN  -  M)P-  EP  +  IQ,  (5.3b) 

=  EP  -  IQ  -  SQ.  (5.3c) 


dN_ 

df 

dP 

dT 

<fQ 

dT 
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The  term  SQ  accounts  for  mortality  during  dispersal.  They  found  that  this  pool  of 
dispersers  was  always  stabilizing,  and  that  the  stabilizing  effect  was  also  produced 
by  a  pool  of  dispersing  prey.  Holt  (1984)  and  Weisser  et  al.  (1997)  extended  these 
results  to  a  system  of  multiple  patches  coupled  through  such  a  dispersal  pool. 

Model  (5.3)  captures  the  essential  fact  that  some  fraction  of  the  predator  popu¬ 
lation  is  dispersing,  and  therefore  not.  consuming  prey  in  habitat  patches.  However, 
like  the  linearly  coupled  Lotka-Volterra.  model  (5,2).  model  (5.3)  makes  some  peculiar 
assumptions  about  the  way  dispersal  occurs.  In  effect,  it  implies  that  there  is  an 
exponential  distribution  of  trip  durations.  Thus  there  is  no  minimum  travel  time,  no 
maximum  travel  time,  and  the  peak  of  the  travel-time  distribution  is  at  zero.  Indeed, 
no  matter  how  long  the  trip,  there  is  a  finite  probability  that  a  given  predator  will 
survive  an  even  longer  trip,  dispersing  without  sustenance. 

Although  the  properties  of  the  dispersal  process  described  by  model  (5.3)  are 
unrealistic,  they  are  no  more  unrealistic  than  other  assumptions  imbedded  in  the 
Lotka-Volterra  model.  Nevertheless,  it  is  important  to  see  if  the  above  stabilizing 
effects  discovered  by  Holt  and  Weisser  et  al.  hold  when  dispersal  is  described  more? 
realistically. 

In  the  next  section,  we  formulate  a  model  similar  to  (5.3)  that  allows  for  an 
arbitrary  distribution  of  trip  durations.  We  show  t  hat  when  only  one  species  disperses, 
tiie  equilibrium  is  almost  always  stabilized  by  including  a  finite  travel  time.  The 
exception  occurs  when  every  trip  lias  exactly  the  same  duration.  In  this  case  there  is 
a  set  of  parameters  values  with  zero- measure  for  which  it  is  not  possible  to  determine 
stability  via  the  linearization  method  we  use.  In  Section  3,  we  formulate  a  two  patch 
version  similar  to  model  (5.2),  and  derive  similar  results.  In  Section  4  we  consider 
multiple  patches  with  two  connection  configurations.  We  have  relegated  some  of  the 
technical  mathematics  required  to  prove  our  results  to  the  Appendix.  We  conclude 
with  a  brief  discussion. 


5.2  Dispersal  Delays  in  1  Patch  Models 

Because  of  the  differences  between  individuals  and  the  vagaries  of  travel,  it  is  rea¬ 
sonable  to  assume  that  dispersal  time  varies  among  individuals  and  between  trips  for 
a  single  individual.  To  incorporate  this  variability,  we  define  a  probability  density 
function,  G{S )  >  0.  for  the  time  it  takes  an  individual  to  disperse,  given  that  the  in- 
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dividual  survives  the  trip,  The  product  C{S)dS  is  the  probability  that  a  successfully 
dispersing  individual  departing  at  time  T  completes  its  trip  between  time  T  +  S  and 
time  T  +  S  +  dS.  Because  each  such  disperser  has  a  nonnegative  travel  time, 

G{S)dS=  1.  (5.4) 


If  there  is  a  constant  probability  per  unit,  time  (M,t)  for  the  disperser  to  perish  while 
travelling,  then  exp (—MaS)  is  the  probability  of  surviving  a  trip  of  duration  S, 

Incorporating  a  distribution  of  travel  times  in  a  single-patch  model  where  both 
prey  and  predators  disperse  gives: 


(IN_ 

~df 

<ip_ 

dT 


( R-AP)N  +  Dn 
[BN  —  M)P+  Dp 


[l 


Gn(S)  e~MfjS  N(T  -  S )  dS  -  N 

r  GP[S)  e~Mps  P[T  -  S)  dS  -  P 

Jo 


(5.5a) 

(5.5b) 


Here,  and  below,  when  the  time  dependence  of  a  variable  is  not  explicitly  indicated 
we  follow  the  convention  that  the  variable  is  evaluated  at  the  current  (undelayed) 
time,  We  assume  that  the  parameters  D Dp,  AIs  and  A/p  are  nonuegative. 

The  analysis  of  model  (5.5a)  is  simplified  by  resealing  variables  and  parameters 
via 


t  =  RT.  $  —  RS,  ft  —  M/R ,  ft.n  =  Mn/R,  ftp  —  Alp/R ,  (5.6a) 

p  =  AP/R ,  n  =  BN/ R,  d.n  =  DN/R ,  dv  =  DP/R.  (5.6b) 


Using  these  new  variables  converts  model  (5.5a)  to  the  dimensionless  form 


h 


P 


(1  -  p)  »  +  dn 
(n  -  ft)  p  +  dv 


s)  ds  —  77. 

ft)  ds  —  p 


(5.7a) 


(5.7b) 


where  fi„{s)  and  gp[s)  arc  the  rescaled  versions  of  Gs[S)  and  GP[S).  The  dot  is  used 
to  denote  a  derivative  with  respect  to  /.  For  the  basic  theory  of  delay  differential 
equations  that  applies  to  model  (5.7)  see  Cushing  (11)77)  and  Kuang  (1993). 

When  both  species  are  mobile,  with  their  own  characteristic  emigration  rate, 
travel-time  distribution,  and  mortality  rate  during  transit,  the  analysis  of  model  (5.7) 
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is  difficult.  For  simplicity,  we  therefore  restrict  attention  to  the  special  cases  in  which 
only  one  species  disperses. 

5.2.1  Predator  Dispersal  ( dn  —  0,  dp  >0) 

When  prey  do  not  disperse  ( dn  —  0),  model  (5.7)  has  two  equilibria.  The  first,  at 
(0,0),  is  always  unstable  to  prey  invasion.  The  second  equilibrium  is  at 

n*  =  H  +  dp(l  -  9p(f*p)h  P*=1-  (5-8) 

where  gp  is  the  (one-sided)  Laplace  transform  of  the  travel-time  distribution  gp.  That 
is,  ^ 

gp(x)=  f  gp{s)  e~xsds.  (5.0) 

JO 

For  real  x,  gp(x)  is  a  positive,  decreasing  function  with  gp( 0)  =  1 . 

We  now  show  that  the  equilibrium  point  (5.8)  is  locally  asymptotically  stable  for 
any  finite  travel-time  distribution  gp{s)  that  lias  measurable  support.  We  begin  by 
linearizing  model  (5.5a)  in  the  neighborhood  of  the  equilibrium  point  (5,8).  Let.  u(t) 
and  c(t)  be  small  perturbations  to  the  equilibrium  point.  That  is,  let 

n(t)  =  n*  +  u(t)'  p(t)=p*  +  v(t),  (5.10) 

with  |u|  <<  n*  and  |u|  <<  p*.  The  dynamics  of  u  and  v  arc  approximately  given  by 
the  linear  system 

it  =  — n*v , 

v  =  p*u  +  dp  |  j  e_fipS3P ( s)  s)  ds  -  gv {pp ) v 

to  which  wc  look  for  exponential  solutions  of  the  form 


Using  (5.12)  in  system  (5.11)  we  obtain  the  system  of  equations  Jw  —  ().  where 

A  n* 

J  =  _  ,  .  (! 

p  A  ~F  dp  [ijp  ( /  tp )  —  gp(pp  4-  A)] 


(5.11a) 

(5.11b) 


(5.12) 
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The  existence  of  a  nontrivial  solution  w  requires  that  det(J) 
the  characteristic  equation 

—  0,  which  in  turn  gives 

11 

(5.14) 

where 

W(A)  =  A[A  +  <«„(/<,,)]  +  nV, 

(5.15) 

and 

A'(A)  =  A  dpgp(fip  +  A). 

(5.16) 

We  next  show  that  fill  roots  of  the  characteristic  e(| nation  (5.14)  have  negative 
real  parts,  and  hence  that  the  equilibrium  (5.8)  is  locally  asymptotically  stable.  To 
do  so  we  first  eliminate  the  possibility  of  roots  with  positive  real  parts  by  assuming 
the  existence  of  such  roots  and  deriving  a  contradiction.  Secondly,  we  eliminate  the 
possibility  of  purely  imaginary  roots,  again  by  deriving  a  contradiction,  Finally,  note 
that  A  —  0  is  not  a  root,  since  //(())  >  0  and  A'(0)  =  0.  The  only  possibility  that 
then  remains  is  that  the  real  parts  of  all  of  t he  roots  are  negative. 

Let  A  =  x  +  iy.  with  x  and  y  real.  Assuming  that  x  >  0,  the  characteristic 
equation  gives 


|//(A)|2  =  |A(A)|2,  (5.17a) 

=  |  {x  +  iy)dpgp{pp  + x +  iy)\2,  (5.17b) 

<  \{x  +  iy)dpgp{ft.p)\2}  (5.17c) 

which,  after  a  little  algebra,  reduces  to 

x4  +  2x2{n*p*  +  if)  +  2xdp(x 2  +  n*p*  +  y2)yp{yp)  +  (if  -  n*p*)2  <  0.  (5.18) 


Since  we  liave  assumed  x  >  0.  each  term  on  the  left-hand  side  of  expression  (5.18) 
is  nonnegativc  and  at  least  one  term  is  positive,  thus  the  left-hand  side  is  posit  ive, 
violating  the  inequality.  As  a  result,  the  roots  of  the  characteristic  equation  cannot 
have  positive  real  parts,  tints  ;r  <  (). 

Now  assume  .r  —  0.  Setting  x  =  0  in  the  characteristic  equation  (5.14),  and 
separating  the  equation  into  its  real  and  imaginary  parts,  shows  that  y  must  be  a 


130 


solution  to  the  system 


n*p*  ~  y 2 

dP  y  9p(Pp) 


dpy 

fOO 

/  gp(s)  e~lip*  sin (ys)  ds , 

Jo 

(5.19a) 

dpy 

r  -x: 

/  9p{$)  e“"pS  cos  {ys)  ds. 

Jo 

(5.19b) 

But  in  the  Appendix,  we  prove  that  for  travel-time  distributions  gp(s)  that  are  finite 
and  have  support  on  a  measurable  set,  there  is  no  real  solution  to  system  (5.19).  Thus 
x  ^  0,  and  since  x  <  0,  it  must  be  that  x  <  0,  and  hence  the  equilibrium  is  locally 
asymptotically  stable. 


Discrete  Delays 

For  our  proof  that  system  (5.19)  has  no  real  solution,  we  must  assume  that  the 
travel-time  distribution  is  finite.  This  requirement  is  satisfied  by  most  biologically 
reasonable  distributions,  including  t  he  exponential  distribution  implicitly  assumed  by 
models  that  use  a  pool  of  dispersers  (Holt  1984,  Weisser  and  Hassell  1996,  Weisser  et 
al.  1997).  or  the  gamma  distribution  that  is  often  used  as  a  convenient  distribution 
because  it  makes  numerical  simulation  easy  (MacDonald  1989). 

However,  if  every  trip  of  every  individual  is  of  exactly  the  same  duration  r,  the 
travel-time  distribution  is  a  delta  function:  yp(s)  =  6(s  —  r).  In  this  case,  the  dis¬ 
tribution  is  not  finite  and  our  proof  does  not  apply.  The  assumption  of  identical 
trips  is  certainly  unrealistic.  Nevertheless,  we  will  analyze  this  ease  below  because 
there  is  a  long  tradition  of  using  discrete  delays  in  population  biology  to  account  for 
individual  development  (Hutchinson  1948;  Wangersky  and  Cunningham  1956.  1957a, 
1957b;  Caswell  1972;  and  many  others),  and  we  would  like  to  compare  the  results  of 
discrete-delays  in  our  model  with  these  results.  Furthermore,  the  analysis  sheds  light 
on  the  reasons  why  dispersal  delays  are  stablizing. 

We  now  show  that,  except  on  a  parameter  set  of  measure  zero,  the  equilibrium 
point  remains  locally  asymptotically  stable  when  all  trips  are  of  the  same  duration. 
The  argument  we  have  already  laid  out,  up  to  and  including  the  characteristic  equa¬ 
tion  (5.14)  holds  with  gp{gv)  =  exp (-//pt).  In  addition,  our  proof  that  there  are  no 
roots  with  positive  real  part  carries  over  to  this  ease.  For  imaginary  roots  A  =  iy, 
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equations  (5. JO)  simplify  to 


(5.20a) 

(5.20b) 


n V  -  V2  =  dpye-w  sin(jrr), 
dpye~(i,,T  =  ( ivye~>lpT  cos(t/t). 


By  (5.20a),  y  ^  0.  Equation  (5.201))  is  satisfied  only  when  cos(yr)  =  1.  in  which  ease 
equation  (5.20a)  gives  y  —  ±y/n*p*  and  equation  (5.20b)  gives, 


T  =  Tk  = 


2kir 


,  for  k  =  0,1,2, 


with  n*  and  p*  given  by  (5.8).  Note  that,  in  (5.8),  n*  depends  upon  r  so  that 
and  (5.21)  must  be  solved  simultaneously  to  find  for  any  set  of  parameters. 


(5.21) 

(5.8) 


When  r  /  r*,.  the  equilibrium  is  locally  asymptotically  stable.  If  r  =  t*.,  then 
the  characteristic  equation  has  purely  imaginary  roots  at  ±iy/v*p*.  For  k  —  0.  the 
model  reduces  to  a  dimensionless  form  of  the  Lot ka- Volt erra  equations  (5.1),  with 
a  neutrally  stable  equilibrium  surrounded  by  a  family  of  periodic  solutions.  When 
k  >  0  we  cannot  infer  the  stability  (or  instability)  of  the  equilibrium  point  from  the 
linearized  analysis. 


5.2.2  Prey  Dispersal  (dn  >  0,  dp  =  0) 

When  only  prey  disperse,  dv  =  0,  and  the  nontrivial  equilibrium  of  model  (5.7) 
I  Jeromes 

n*  =  p,  P*  =  1  -  dn(l  -  finiUn))-  (5-22) 

If  either  the  prey  emigration  rate  ((/„)  or  the  prey  mortality  rate  while  dispersing  (//.„) 
is  too  large,  both  prey  and  predators  are  unable  to  persist,  and  the  equilibrium  at 
(0,  U)  becomes  stable. 

Linearizing  about  (n*,p*),  with  p*  assumed  positive,  gives  the  same  characteristic 
equation  (5.14)  as  in  the  predator  dispersal  case,  with  dp,  pp,  and  gp  replaced  by 
</„,  and  gn.  Thus  th<'  results  of  Section  2.1  hold  for  mobile  prey  and  sedentary 
predators  whenever  p*  >  0. 


132 


5.3  Dispersal  Delays  in  2  Patch  Models 


Incorporating  a  distribution  of  travel  times  for  both  the  predator  and  its  prey  in  the 
2- patch  Lotka-Volterra  model  {5.2}  gives 

dNi/dT  =  (/?  -  APi)  Ni  +  Dn  |jf  G*{S)  e“M*s  Nj{T  -  S)  dS  -  N,  , 

(5.23a) 

dPi/dT  =  {BNi-  M)  P  +  DP  [jf  GV(5)  e“MpS  P^T  -S)dS-  Pt  , 

(5.23b) 


for  i,  j  =  1, 2  and  j  ^  *.  Using  the  rescaled  variables  (5.6).  now  subscripted,  converts 
model  (5.23)  to  the  dimensionless  form 


n,  =  (1  -  pj)  ttj  +  ^  pn(s)  e  nj(t  -  s)  ds  -  rq  , 

Pi  =  -  p.)  pi  +  dp  flp(a)  e"^ Pj{/  -  s)  ds  -  pi  . 


(5.24a) 

(5.24b) 


Again,  we  analyze  only  the  cases  in  which  one  species  disperses.  When  only  preda¬ 
tors  disperse  (i.e.,  d„  =  0,  dp  >  ()},  model  (5,24)  has  a  unique  positive  equilibrium 
that  is  spatially  homogeneous  with  densities  equal  to  the  equilibrium  densities  in  the 
1 -patch  model  with  predator  dispersal:  ri\  =  n‘2  =  n*  and  p*  =  p$  =  p*  with  n*  and 
p*  given  by  (5.8).  Linearizing  model  (5.24)  in  the  neighborhood  of  this  equilibrium 
gives,  in  analogy  with  (5.13) 


A  B 
B  A 


(5.25) 


wit  h 


A  = 


A  n* 

—p  A  +  dpQpipp) 


and  B  - 


0  0 
0  dp(jp(fip  -(-  A) 


(5.26) 


Setting  det(J)  =  0  yields  the  following  characteristic  equation: 


H2(X)  =  A’a(A), 


(5.27) 


with  H{ A)  and  A'(A)  as  in  the  single  patch  model  (i.e.,  given  by  (5.15)  and  (5.16)). 


The  arguments  in  Section  2  can  now  be  applied  almost  unchanged.  For  any 
finite  travel-time  distribution  with  support  on  a  measurable  set,  the  equilibrium  is 
stabilized.  The  only  change  comes  in  the  discrete-delay  case,  when  the  travel-time 
distribution  is  a  delta  function.  In  this  case  equations  (5.20)  must  be  modified, 
because  any  A  —  iy  that  satisfies  H{ A)  =  A'(A)  or  H( A)  =  —A' (A)  is  a  purely 
imaginary  eigenvalue.  Thus  equations  (5.20)  become 

n*p*  —  y2  =  ±dp  y  .sin(yr).  (5.28a) 

dp y e-#tpT  -  ±dpye~tlp7  cos(yr),  (5.28b) 


which  now  give  y  =  and  r  =  77-/2  (cf.  equation  (5.21)).  When  r  is  equal  to 

one  of  these  special  values,  the  eigenvalues  are  purely  imaginary,  and  thus  linearization 
cannot  be  used  to  determine  the  stability  of  the  equilibrium  point. 

When  only  prey  disperse  (i.e,,  <l,x  >  0,  dp  =  0),  the  unique  nont  rivial  equilibrium  is 
also  spatially  homogeneous  and  equal  to  the  equilibrium  densities  of  the  1- patch,  prey- 
dispersal  model  {equations  (5.22)).  which  we  assume  to  be  positive.  The  linearization 
still  gives  the  matrix  (5.25),  but  with  B  now  given  by 


B  = 


— d,t(}j\ Uhi  +  A)  0 
0  0 


(5.2SJ) 


Setting  det(J)  —  0  again  gives  the  characteristic  equation  (5.27).  but  with  all  sub¬ 
scripts  changed  from  p  to  n.  All  of  the  just  derived  two- patch  predator-dispersal 
results  therefore  carry  over  to  the  prey-dispersal  case. 


5.4  Multiple  Patches 


We  are  also  able  to  prove  that  one-species  dispersal  delays  are  stabilizing  in  con¬ 
figurations  of  an  arbitrary  number  of  patches  that  admit  a  spatially  homogeneous 
equilibrium.  In  such  configurations  every  patch  is  identical  to  every  other  patch.  A 
ring  of  m  >  3  patches  is  perhaps  the  simplest  example,  hi  the  predator  dispersal 
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case,  the  ring  model  is  given  by: 


hi  =  (1  —  Pi)  rii 

Pi  =  (ni-ii)pt  + 

Y  9P(s)  e~^s  (pi_i(t  -  s)  +  pi+i  ((  -  a))  ds  -  2 p,  , 


(5.30a) 

(5.30b) 


for  i  =  1, . . . ,  m.  To  close  the  ring,  define  po  =  pm  and  p„l+  j  =  p\ . 


The  equilibrium  (5.8)  remains  unchanged.  Linearizing  around  it  and  substituting 
an  exponential  solution  gives  the  characteristic  equation  det(J{m))  =  0,  with  the 
2m  x  2m  matrix  J{m)  given  by 


J(m)  = 


A 

0 

0 

& 

A 

0 

0 

0 

& 

A 

0 

0 

0 

& 

A 

0 

. . . 

0 

A 

(5.31) 


with  blocks  A  and  |B  given  by  (5.26). 

J(m)  is  an  example  of  a  block-circulant  matrix  with  2x2  blocks.  That  is,  it  has 
the  form 

Ad  A[  •••  Ar„_i 

C  =  I  Am"'  A .  A"-*  ,  (5.32) 

A i  A •>  A(| 

where  the  A;  are  2x2  matrices.  The  determinant  of  such  matrices  is  given  by 

1 


det  C=I]  detTf, 


(5.33) 


f=0 


with 


m— 1 


j=0 


(5.34) 


(see,  for  example,  Friedman  (1961,  Theorem  6)).  Applying  these  formulae  to  Jtm) 


gives 


If  X  is  an  eigenvalue,  at  least  one  term  in  the  product  (5.35)  will  vanish,  i.e., 

H{ A)  =  cos  K(X)  (5.36) 

for  some  f  <  m.  But  taking  absolute  values  of  both  sides  and  squaring  gives  |//(A)|*  < 
|A'(A)|2.  which  brings  ns  back  to  equal  ion  (5.17)  and  the  single  patch  case.  Distributed 
dispersal  delays  are  always  stabilizing  in  this  “ring"  model. 

For  discrete  delays  of  duration  r.  the  real  and  imaginary  parts  of  equation  (5.36) 
give,  in  analogy  to  equations  (5.20). 


While  y  =  0  is  always  a  solution  to  (5,37b),  it  is  never  a  solution  to  (5.37a).  We 
therefore  take  ,</  ^  0.  If  the  number  of  patches  in  the  ring  (m)  is  odd,  there  is  only 
one  value  of  f  with  f  <  m  for  which  (5.37b)  has  a  solution:  it  is  f  =  0.  In  this  case 
purely  imaginary  eigenvalues  occur  at  r  =  r*.,  with  given  by  equation  (5.21).  If 
the  number  of  patches  is  even,  imaginary  eigenvalues  occur  at  r  -  r*/2.  For  k  even 
they  arise  from  (5.37b)  with  t  =  0.  For  k  odd,  they  come  from  (5.37b)  with  (  =  m/2. 
Linearization  is  uninformative  for  these  delays. 

Another  configuration  with  all  patches  identical  is  obtained  if  each  patch  is  coupled 
to  every  other  patch  in  exactly  the  same  way  (e.g.,  through  a  pool  of  dispersers).  The 
resulting  model  in  this  case  is 


n, 

Pi 


(1  -  pi)  Hi 

(«i  “  k)  Vi  + 


(5.38a) 


m  -  1 

■OC 


X]  f  f  9] >(s)  e  Wpjit  -  s )  ds  -  (m  -  l)pf 1 ,  (5 

j-u&i Uo  J  J 


.38b) 


for  i  =  1 . m.  For  m  =  2,  this  model  reduces  to  the  two-patch  model  of  Section 

2.1:  for  in  =  3  it  is  the  same  as  model  (5.30).  We  therefore  take  m  >  -1. 
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Linearization  around  the  equilibrium  {again  given  by  (5.8))  gives 


(5.39) 


A 


and  the  resulting  characteristic  equation  is 


det(J,m))  =  []  H{ A)  - 


777  -  1 


1 


K{\)  .  (5.40) 


By  the  argument  used  on  the  determinant  (5.35),  it  follows  that,  for  a  distributed 
delay,  the  conclusions  are  the  same  as  in  the  single-patch  case.  A  discrete-delay  of 
duration  r  is  also  stabilizing,  unless  r  =  r*. 

5.5  Discussion 

The  predator-prey  models  we  formulated  above,  wherein  one  species  disperses  between 
habitat  patches  while  t  he  other  does  not,  all  show  that  a  dispersal  delay  almost  always 
stabilizes  the  spatially  homogeneous  positive  equilibrium.  It  is  well  known  that,  the 
inclusion  of  a  delay  can  lead  to  a  qualitative  change  in  the  dynamics  of  a  model,  but 
it  is  typically  the  case  that  an  increase  in  the  delay  produces  instability  and  gives  rise 
to  stable  periodic  solutions  (see,  e.  g.,  MacDonald  1989,  p.  8,  15).  In  this  sense,  our 
results  can  be  seen  as  counterintuitive,  although  some  models  with  delay-dependent 
parameters  exhibit  stability  switches  from  stable  to  unstable  and  back  to  stable  again 
as  the  delay  increases  (see,  e.  g.  Beretta  and  Kuang  2001). 

What  is  the  mechanism  by  which  dispersal  delays  act  to  stabilize  these  systems? 
First  let  us  say  that  certain  mechanisms  are  not  responsible.  The  stabilizing  effect 
is  not  (strictly  speaking)  a  metapopulation  effect.  After  all.  it  is  evident  in  the  one- 
patch  model.  It  is  not  an  effect  of  a  cost  of  dispersal  (in  terms  of  increased  mortality). 
After  all.  the  mechanism  operates  when  ftp  or  ft,,  vanish.  Finally  it  is  not  a  result  of 
spatial  heterogeneity;  every  patch  is  identical  to  every  other  patch  in  our  models. 

The  stabilizing  mechanism  that  does  operate  is  evident  in  the  one-patch  predator- 
dispersal  model.  In  particular,  consider  the  case  where  /ip  =  0.  In  this  case,  the  term 
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that  accounts  for  predator  dispersal  can  be  rewritten  as 


p(0 


'irfip(s)p[t-s)ds 

P(t) 


(5.41} 


The  term  between  the  curly  brackets  represents  the  instantaneous  per  capita  net  mi¬ 
gration  rate  (i.e.,  immigration  minus  emigration).  It  is  density  dependent.  If  the 
population  in  the  patch  is  high  relative  to  historically  average  values  (the  average 
being  taken  with  respect  to  the  weighting  function  gp(s))>  then  net  migration  is  neg¬ 
ative:  if  the  current,  population  size  is  low  relative  to  historical  averages,  then  the 
net  migration  rate  is  positive.  This  term  therefore  has  t  he  effect  of  damping  oscilla¬ 
tions  and  enhancing  stability.  The  same  mechanism  also  works  in  the  multiple-patch 
scenario. 


Murdoch  et  al.  (1992)  also  showed  that  “temporal  density-dependence”  in  immi¬ 
gration  rates  can  stabilize  predator-prey  metapopulation  dynamics.  The  fundamental 
difference  between  their  results  and  ours  is  the  mechanism  by  which  the  density- 
dependence  is  generated.  In  their  model,  spatial  heterogeneity  in  the  demographic 
parameters  generates  asynchronous  population  dynamics  between  connected  habitat 
patches,  which  leads  to  a  decoupling  of  local  immigration  rates  from  local  popula¬ 
tion  density.  In  our  model,  no  such  spatial  heterogeneity  is  required  to  generate  the 
decoupling.  Instead,  dispersal  tends  to  synchronize  the  dynamics  between  patches, 
while  the  delay  decouples  immigration  rates  from  local  density. 

When  the  delay  is  discrete,  there  is  a  set  of  delays  for  which  we  have  not  been 
able  to  determine  the  stability  of  the  equilibrium.  Because  this  set  of  delays  lias  zero 
measure,  these  cases  are  biologically  irrelevant.  Nevert  heless,  an  understanding  of  the 
dynamics  in  these  cases  would  complete  the  mathematical  analysis.  We  conject  ure 
that  for  these  special  values  of  t  he  delay,  dispersal  does  not  stabilize  the  equilibrium. 

The  models  we  have  analyzed,  while  more  complex  than  the  Lot  ka-Vol terra  model, 
are  still  simple  in  the  extreme.  An  important  simplification  we  have  made  is  that  ev¬ 
ery  patch  is  identical  to  every  other  patch  including  being  connected  via  dispersal  to 
the  same  number  of  equidistant  patches.  In  real  systems,  this  assumption  is  violated. 
When  the  distance  between  patches  is  not  constant,  the  travel  time  distribution  will 
differ  for  different  pairs  of  patches.  Models  of  this  type  are  notoriously  difficult  to 
analyze,  but  it  is  important  to  know  the  extent  to  which  our  results  rely  on  this  as¬ 
sumption.  The  interaction  between  more  realistic  population  dynamics  and  dispersal 
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delays  is  a  topic  that  we  are  currently  investigating  and  will  report  elsewhere. 
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Abstract 


It.  takes  time  for  individuals  to  move  from  place  to  place.  This  travel  time  can 
be  incorporated  into  metapopulation  models  via  a  delay  in  the  interpatch  migration 
term.  Such  a  term  has  been  shown  to  stabilize  the  positive  equilibrium  of  the  classical 
Lotka-Volterra  predator  prey  system  with  one  species  (either  the  predator  or  the 
prey)  dispersing. 

We  study  a  more  realistic,  Rosen  zweig- Mac  Arthur,  model  that  includes  a  carrying 
capacity  for  the  prey,  and  saturating  functional  response  for  the  predator.  We  show 
that  dispersal  delays  can  stabilize  the  predator  prey  equilibrium  point  despite  the 
presence  of  a  Type  II  functional  response  that  is  known  to  be  destabilizing.  We  also 
show  that  dispersal  delays  reduce  the  amplitude  of  oscillat  ions  when  the  equilibrium 
is  unstable,  and  therefore  may  help  resolve  the  paradox  of  enrichment.. 


6.1  Introduction 

The  basic  models  of  predator  prey  and  host  parasitoid  systems  predict  unstable  equi¬ 
libria.  often  accompanied  by  large-amplitude  oscillations  in  both  species.  These  oscil¬ 
lations  drive  the  populations  to  low  densities,  and  have  been  interpreted  as  potential 
causes  of  extinction.  In  contrast,  natural  predator  prey  systems  seem  to  persist  for 
long  periods.  Theoreticians  and  experimentalists  have  suggested  a  number  of  p<e 
tential  processes  that  might  resolve  this  conflict  between  models  and  data  (sec.  for 
example,  May,  1973:  Hassell,  1978;  Crawley.  1992;  Mueller  C  Joshi,  2000).  Spatial 
processes,  and  in  particular  metapopulation  structure,  have  garnered  significant  at¬ 
tention  (Taylor,  1990:  Briggs  &  Hoopes,  2004). 

Dispersal,  the  process  that  distinguishes  spatial  models  from  their  nonspat  ial  coun¬ 
terparts,  has  been  added  to  predator-  prey  models  in  many  different  ways,  with  vary¬ 
ing  effects  on  stability  (Briggs  &  Hoopes,  2004).  One  way  to  include  dispersal  is  to 
distinguish  a  class  of  dispersing  individuals,  that,  while  dispersing,  do  not  participate 
in  tb<’  predator-prey  interaction.  A  number  of  authors  have  shown  that  including 
such  a  pool  of  dispersers  (be  they  predators  or  prey)  in  a  Lotka  Volterra  model  st  a¬ 
bilizes  coexistence  at  an  equilibrium  point.  The  models  of  Holt  (1984),  Weisser 
Hassell  (1996)  and  Weisser  ct  al.  (1997)  include  the  dispersal  pool  explicitly,  and 
couple  it  to  the  dynamics  within  a  patch  via  constant  per  capita  immigration  and 
emigration  rates.  These  models  implicitly  assume  an  exponential  distribution  of  the 
t  ime  that  an  individual  spends  dispersing. 

Exponential  travel-t  ime  distributions,  however,  have  some  biological  peculiarities. 
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For  example,  there  is  no  maximum  travel  time,  and  the  modal  travel-time  is  zero. 
To  see  if  these  implicit  assumptions  play  a  role  in  stabilizing  the  equilibrium,  Neu- 
bert  et  at.  (2002)  relaxed  this  assumption  by  prescribing  an  arbitrary  distribution 
of  dispersal  times.  They  showed  that,  except  in  cases  so  rare  as  to  be  biologically 
irrelevant,  the  stabilizing  effect  of  such  “dispersal  delays”  remains. 


All  of  these  analyses  are  based  upon  the  Lotka  Volterra  predator  prey  model 

dN 

—  =  {R-AP)N,  (6.1a) 

dP 

—  =  (BN-M)P,  (6.1b) 

dT 

where  N  is  the  population  density  of  the  prey  and  P  is  the  population  density  of 
the  predator.  The  prey  population  has  a  constant  per  capita  growth  rate  /?.  and  the 
predator  population  has  a  constant  per  capita  mortality  rate  M.  The  predator-prey 
interaction  is  captured  bv  linear  functional  and  numerical  responses,  scaled  by  the 
parameters  A  and  B.  The  parameters  /?,  A.  B.  and  A/  are  assumed  to  be  positive. 


Model  (6.1)  has  a  unique  coexistence  equilibrium  point  (i.e.,  an  equilibrium  point 
at  which  both  species  have  positive  densities)  at  N  —  M/B ,  P  —  R/A.  This  equi¬ 
librium  point  is  a  center,  surrounded  by  a  family  of  periodic  orbits  whose  amplitudes 
depend  on  the  initial  population  sizes.  Adding  either  predator  or  prey  dispersal  to 
this  model  stabilizes  the  equilibrium  point  if  dispersal  delays  are  accounted  for  (Neu- 
bert.  et  a/.,  2002).  In  the  absence  of  delays,  predator  dispersal  reduces  the  amplitude 
of  the  oscillations  but  does  not  stabilize  the  equilibrium  point  (Jansen,  1905:  Jansen 
&i  de  Roos,  2000).  Increasing  the  number  of  patches  in  this  model  gives  rise  to  other 
equilibria  in  which  the  prey  are  absent  front  one  or  more  patches  (see,  for  example, 
Feng  k  Hinson,  2005)  that  we  do  not  consider  here. 


Model  (6.1).  and  its  spat  ial  extensions,  have  been  criticized  as  being  oversimplified 
for  two  reasons.  First,  in  the  absence  of  the  predators,  the  prey  grow  exponentially 
without  bound.  Second,  the  per  capita  rate  of  consumption  of  prey  by  predators 
grows  in  proportion  to  the  prey  population  size,  implying  that  individual  predators 
can  process  prey  items  infinitely  fast.  These  faults  are  eliminated  in  the  Rosenzweig- 
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MacArthur  model  (Rosenzweig  &  MacArthur,  1963) 


(IN_ 

df 

dP 

dT 


(6.2b) 


for  the  predators  that  results  in  a  saturating  functional  response.  Here,  A  is  the 
maximum  rate  at  which  an  individual  predator  can  consume  prey  and  C  is  the  prey 
density  at  which  an  individual  predator’s  consumption  rate  equals  A/2.  The  ratio 
D/A  gives  the  fraction  of  consumed  prey  that  are  converted  into  predators. 

The  dynamics  of  model  (6.2)  an>  more  complicated  than  those  of  model  (6.1) 
(Kot.  2001).  For  small  values  of  carrying  capacity,  the  coexistence  equilibrium  point 
is  locally  asymptotically  stable.  As  the  carrying  capacity  increases  beyond  some 
threshold  value,  the  equilibrium  point  becomes  unstable,  and  trajectories  are  drawn 
onto  a  single  stable  limit  cycle.  The  amplitude  of  predator-prey  oscillations  increases 
with  increasing  prey’s  carrying  capacity,  reaching  vanishingly  small  densit  ies  at  which 
natural  populations  cannot  persist.  This  destabilization  by  increasing  prey  carrying 
capacity  is  known  as  the  'paradox  of  enrichment '  (Rosenzweig,  1971;  May,  1972; 
Gilpin,  1972). 

Here  we  present  three  major  findings.  First,  we  show  that  dispersal  delays  can  sta¬ 
bilize  the  coexistence  equilibrium  point  of  model  (6.2)  (as  they  did  in  model  (6.1))  by 
delineal  ing  the  stability  region  in  parameter  space.  Second,  we  show  that  for  many  pa¬ 
rameter  values,  stability  persists  in  a  so-called  “Type  II  model"  wherein  prey  growth 
is  density  independent  (i.e.  model  (6.2)  in  the  limit  of  infinite  carrying  capacity  A*). 
We  thus  establish  that  delayed  dispersal  can  overcome  a  destabilizing  Type  H  func¬ 
tional  response  even  in  the  absence  of  stabilizing  prey  density-dependence.  Finally, 
we  show  t  hat  dispersal  delays  help  resolve  the  paradox  of  enrichment  by  reducing  the 
amplitude  of  oscillations  when  the  equilibrium  is  unstable,  thereby  preventing  the 
small  population  sizes  that  might  lead  to  extinction. 

We  begin,  in  the  next  section,  by  constructing  a  Rosenzweig- Mat1  Arthur  mode! 
that  incorporates  dispersal  delays.  Using  the  methods  outlined  in  Neubert  et  ai 
(2002),  it  can  be  shown  that  if  dispersal  delays  stabilize  t  he  single  patch  model  they 
also  stabilize  a  spat  ially  homogeneous  equilibrium  of  a  model  with  an  arbitrary  num¬ 
ber  of  identical  patches.  Therefore,  we  limit  our  investigation  to  a  single  habitat 


patch  from  which  only  predators  disperse.  We  then  present  results  for  two  types  of 
dispersal  delay:  a  discrete  delay  that  implies  that  all  individuals  spend  exactly  the 
same  amount  of  time  away  from  the  patch,  and  a  distributed  delay  that  accounts  for 
differences  in.  for  example,  dispersal  ability  between  individuals.  For  discrete-delays 
our  results  are  derived  from  numerical  simulations.  In  the  case  of  a  distributed  delay 
with  Erlang  distribution,  we  analytically  derive  a  polynomial  characteristic  equation, 
whose  roots  we  find  numerically.  We  conclude  with  a  brief  discussion. 


6.2  Model 


The  model  that  we  analyze, 


dN_ 

dT 

dP 

dT 


ANP 
C  +  N' 


BNP 
C  +  N 


-  NfP  +  D 


G(S)  e~MpsP(T  -S)dS-P 


(6.3a) 

(6,3b) 


describes  the  dynamics  of  a  sedentary  prey  and  a  mobile  predator  in  a  single  habitat 
patch.  Individual  predators  emigrate  from  the  patch  at  the  constant,  per  capita  rate 
D,  and  return  S  units  of  time  after  their  departure.1  To  account  for  the  differences 
in  dispersal  abilit  ies  between  predators,  we  define  a  distribution  of  dispersal  delays, 
G{S)  >  0,  for  the  time  a  predator  takes  to  disperse,  given  that  it  survives  the  trip 
(Neubert,  et  ai,  2002).  Because  all  dispersal  times  are  non  negative  it  follows  that 
f™  G(S)  dS  =  1  (see  also  Azer  k  van  den  Driessche,  2006).  We  assume  that  the 
probability  of  surviving  a  trip  of  duration  S  is  e~MpS,  where  Mp  is  the  mortality  rate 
during  the  migration. 

Model  (6.3)  takes  the  form  of  a  delay  differential  equation  with  distributed  delay. 
For  examples  of  how  such  equations  have  been  used  in  other  types  of  ecological  models, 
and  for  how  they  may  be  analyzed,  see  the  books  by  Kuang  (1993)  and  MacDonald 
(1989). 

In  order  to  reduce  the  number  of  parameters,  and  simplify  our  analyses,  we  scale 


’For  notations]  convenience,  a  variable  with  no  time  dependence  explicitly  given  is  to  be  evaluated 
at  the  current  (undelayed)  time, 


147 


the  variables  and  parameters  of  the  model  (6.3)  according  to 


f  —  RT ,  s  —  RS ,  ft M/R,  d  =  D/R,  pp  —  Mp/R, 
p  =  AP/RC,  n  =  DN/RC ,  £  =  R/D.  k  =  KB/RG. 


Substitution  into  system  (6.3)  gives  the  dimensionless  form 


h  = 
V  = 


where  g{s)  is  the  scaled  version  of  G(S). 


(6.4a) 

(6.4b) 


(6.5a) 

(6.5b) 


In  See.  3.  we  focus  on  the  effects  of  n  (the  dimensionless  carrying  capacity) 
and  d  (the  dimensionless  emigration  rate)  on  the  stability  of  the  unique  coexistence 
equilibrium  for  model  (6.5): 


/<  +  d(l  -  g{fip)) 

1  -£[//  +  d(l  -  ()(//,,))]' 


P*  =  ^1  -  ^7^  (1  +  en*) 


(6.6) 


Here,  g(x)  is  t ht'  (one-sided)  Laplace  transform  of  the  travel-time  distribution  //(*), 


The  equilibrium  (6,6)  is  positive  only  if 


(6.7) 


K  (1  -  Eft)  -  p 

(k£+  1)  [1  -g(nP)]' 


(6.8) 


Ifd  is  too  large,  and  inequality  (6.8)  is  violated,  the  predators  do  not  spend  sufficient 
time  feeding  on  the  prey  patch  to  maintain  a  positive  growth  rate  and  are  extirpated 
as  a  result. 


To  determine  the  stability  of  the  coexistence  equilibrium  point  (6.6)  of  model  (6.5) 
we  must  determine  the  fate  of  small  perturbations,  u(t)  and  v(t)<  to  the  coexistence 
equilibrium,  Sot 

7i  (/)  =  n *  +  u{t),  p  (f)  —  p*  +  v(i).  (6.9) 

For  )«|  and  |u|  sufficiently  small,  the  dynamics  of  these  perturbations  are  approxi- 
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mated  by  the  linear  system 
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Looking  for  solutions  to  (6.10)  of  the  form 


=  w  ^  0, 


we  find  that  A  and  w  must  satisfy 


(J  —  AI)  w  —  0, 


(6.10a) 


(6.10b) 


(6.11) 


(6.12) 


where  J  is  the  .Jacobian  matrix 


.1 
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k  (l+sn*)2 

- V‘  . 

(1+en*)1 
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d[g{fi.p  +  A)  -  g{nP)] 


(0.13) 


Equation  (6.12)  1ms  solutions  with  w  ^  0  only  if  del  (J  —  AI)  =  0,  which  translates 


to 


H{  A)  =  K(  A), 


(6.14) 


with 


H{\) 
K(  A) 
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+  k  +  (1  +  en*)2 
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k  +  (l  +  en*)2 


-  1 


[A  +  dy(/q,)]  + 


w*p* 

(1  +  en*)3 


dgiiip  +  A). 


(6.15a) 

(6.15b) 


The  roots  of  this  “characteristic”  equation  are  the  eigenvalues;  they  are.  in  general, 
complex  numbers. 

The  real  parts  of  the  eigenvalues  determine  the  stability  of  the  equilibrium  point. 
If  all  of  the  eigenvalues  have  negative  real  parts,  u  and  v  will  vanish  in  the  limit 
t  —>  oo,  and  the  equilibrium  point  is  therefore  locally  stable.  If  any  eigenvalue  has  a 
positive  real  part,  the  perturbations  grow,  and  the  equilibrium  is  unstable.  Note  that 
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(6.14)  and  (6.15)  together  imply  that  A  —  0  is  not  an  eigenvalue  since  n*  and  p*  are 
positive. 

hi  the  absence  of  dispersal,  d  =  0.  In  this  case,  local  stability  of  the  equilibrium 
point  is  guaranteed  from  (6.8)  and  (6.14)  if 


1 1  1 +  EU 

— 1 -  <  K  <  — - 

\~ep  e(l-Ep) 


(6.16) 


For  finite  k,  if  the  left-hand  inequality  is  violated  the  predator  is  extirpated,  since  for 
this  parameter  range  there  is  no  positive  steady  state;  see  inequality  (6.8).  Violation 
of  the  right-hand  inequality  results  in  a  Hopf  bifurcation  and  a  predator-prey  limit 
cycle  (Kot,  2001).  Note  that  in  the  limit  k  —>  oc,  the  equilibrium  point  is  never 
stable. 


6.3  Results 

6.3.1  Discrete  travel  time 

If  the  duration  of  every  dispersal  event,  of  every  individual  is  exactly  r.  then  t  he  dis¬ 
persal  delay  distribution  is  a  delta  function:  <](s)  —  S(s— r)  and  </{.r)  =  exp{—  rx).  We 
have  been  unable  to  analyt  ically  infer  the  local  stability  of  the  coexistence  equilibrium 
in  this  case,  as  the  characteristic  equation  {6. 14)  is  a  transcendental  equation  with 
infinitely  many  solutions.  Therefore,  we  illustrate  our  results  (in  Fig.  6-1)  using  nu¬ 
merically  generated  stability  diagrams  in  the  (r,  d)  parameter  plane  for  various  values 
of  h.  For  each  combination  of  the  parameters,  we  (i)  calculated  the  equilibrium  point 
(6.6).  (ii)  for  coexistence  equilibria  we  chose  a  random  initial  condition  for  the  prey 
and  the  predator  uniformly  distributed  between  50%  and  150%  of  the  equilibrium 
values,  (iii)  using  the  Simulink  package  in  Matlab,  we  simulated  the  model  (6.5)  and 
discarded  the  transient  dynamics.  We  then  distinguished  three  sets  in  (r. d)  parame¬ 
ter  plane:  (a)  a  set  of  parameters  for  which  the  coexistence  equilibrium  does  not  exist 
(because  inequality  (6.8)  is  violated),  (b)  a  set  for  which  the  coexistence  equilibrium 
exists  but  it  is  unstable,  and  (e)  a  set  for  which  the  coexistence  equilibrium  is  stable. 

We  start  (Figure  6- 1 A)  with  a  ease  that  is  stable  in  the  absence  of  the  disper¬ 
sal  delay  (i.e.,  satisfying  (6.16)).  As  expected,  the  coexistence  equilibrium  is  stable 
everywhere  it  exists.  In  Figures  6-1B-D  the  values  of  k  violate  the  right-most  in¬ 
equality  in  (6.16).  For  these  values  of  the  carrying  capacity,  the  equilibrium  point  of 
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A 


travel  time  (t) 


D 


travel  time  (t) 


Figure  6-1:  Stability  diagrams  from  simulations  of  model  (6.5)  with  discrete  dispersal 
delay  for  various  values  of  k,  e  =  0.01,  /r  =  ftp  =  1:  A)  k  ~  50;  B}  k  =  150; 
C)  k  =  500;  D)  k  —  5000.  White  areas  designate  a  stable  equilibrium  point,  dark 
gray  stands  for  an  unstable  equilibrium  point  and  an  area  where  t  here  is  no  positive 
coexistence  equilibrium  is  shown  in  light  gray,  bounded  by  the  black  curve,  i.e.,  the 
ease  of  equality  (6.8). 
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the  Rosenzweig-MacArthur  model  (6.2)  is  unstable,  surrounded  by  a  predator- prey 
cycle.  Dispersal  delay  dampens  the  predator-prey  oscillation  resulting  in  the  area  of 
stable  equilibrium  shown  in  white.  For  k  1000,  the  stability  region  reduces  to  four 
“islands" . 

Since  stability  diagrams  do  not  reveal  the  details  of  the  unstable  behavior,  we 
generated  bifurcation  diagrams  for  different  values  of  parameters,  for  both  predator 
and  prey  densities.  All  diagrams  exhibit  qualitatively  similar  behavior,  so  we  show 
only  one  bifurcation  diagram  for  prey  density  with  r  as  the  bifurcation  parameter 
(Fig.  6-2). 

For  each  value  of  r  we  simulated  model  (6.5),  discarded  the  transient  dynamics, 
and  present  only  the  final  behavior  by  plotting  only  the  local  maxima  and  minima  of 
the  trajectory.  Stable  equilibria  therefore  appear  as  a  single  point..  Oscillations  wit  h 
one  peak  appear  as  two  points,  and  oscillations  with  two  peaks  appear  as  four  points, 
et  cetera.  Quasi-periodie  and  aperiodic  oscillations  appear  as  “smears." 

In  addition  to  quasi-periodie  and  aperiodic  behavior,  the  bifurcation  diagrams 
also  reveal  the  coexistence  of  multiple  attractors.  In  Fig.  6-2A.  we  increased  r  from 
0  to  7  in  small  steps,  using  the  end  of  the  simulation  for  one  value  of  r  as  the 
initial  condition  of  the  simulation  for  the  following  value  of  r.  We  followed  the  same 
procedure  in  Fig.  6-2B.  except  that  we  decreased  r  from  7  to  0.  For  values  of  r  in 
the  shaded  regions  of  Fig.  6-2.  solutions  converge  to  different,  attractors  depending 
on  initial  conditions. 

In  Fig.  6-6  we  categorize  the  dynamics  of  the  Type  II  model, 
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(6.17a) 

(6.17b) 


in  the  (r.  d)  parameter  space  over  a  range  of  e.  In  Fig.  6-3A,  the  predat  or's  funct  ional 
response  is  strong  (e  is  relatively  large)  and  the  equilibrium  cannot  be  stabilized  by 
dispersal  delays.  As  e  decreases,  however,  stable  islands  grow  in  number  and  in  size. 
In  the  limit  e  — *  0,  Ncubcrt  et  al.  (2002)  showed  the  equilibrium  is  stable  everywhere 
except  for  a  set  of  measure  zero  in  t  he  (r,  d)  plane.  Comparing  Fig,  6-3C  with  Fig.  6- 
1D  shows  that  the  stability  properties  of  the  Type  II  model  are  essentially  the  same 
as  the  Rosenzweig-MacArthur  model  with  large  carrying  capacity. 

The  stability  in  the  Type  II  model  implies  that  dispersal  delays  can  help  resolve 
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Figure  (i-2:  Bifurcation  diagrams  for  the  prey  population  density  of  model  (6.5)  with 
discrete  dispersal  delay  (minimum  and  maximum  population  densities),  d  =  35.  e  = 
0.01,  k  =  5000,  ft.  =  ftp  =1.  A)  r  is  changed  forwards;  B)  r  is  changed  backwards. 
Shaded  regions  depict  the  coexistence  of  multiple  attractors.  Detailed  explanation  in 
the  text. 
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Figure  0-3:  Stability  area  of  the  Type  II  model  (6.17)  with  discrete  dispersal  delay  for 
various  e:  A)  e  =  U.2:  B)  e  =  0.1;  C)  e  =  0.01;  and  D)  e  =  0.0002.  Other  parameter 
values  are  as  in  Fig.  1. 


Figure  6-4:  Bifurcations  diagram  for  model  (6.6)  with  k  as  bifurcation  parameter  show 
that  the  amplitude  of  oscillation  is  significantly  smaller  in  the  presence  of  dispersal 
delays;  d  =  35.  e  —  0.01,  ft  =  fip  =1.  A)  prey  population,  no  dispersal  (r  =  0);  B) 
prey  population,  t  —  3:  C)  predator  population,  no  dispersal  (r  =  0);  D)  predator 
population,  r  =  3. 


the  paradox  of  enrichment.  In  the  MacArthur-Rosenzweig  model  without  dispersal, 
the  amplitude  of  oscillation  increases  with  increasing  carrying  capacity  and  the  pop¬ 
ulation  soon  reaches  vanishingly  small  densities.  We  illustrate  this  with  a  bifurcation 
diagram  with  k  as  a  bifurcation  parameter  (Fig.  6-4 A,  C).  We  again  present  only  the 
long-term  dynamics  by  plotting  local  minima  and  maxima  of  the  trajectory.  The  min¬ 
imal  population  density  decreases  rapidly  with  increasing  k  and  eventually  becomes 
dominated  by  numerical  round-off  errors,  so  the  graphs  in  Fig.  6-4A  and  C  appear 
blurred.  For  comparison,  in  Fig.  6-4  B  and  D  we  show  how  the  amplitude  of  the  os¬ 
cillation  changes  with  increasing  capacity  in  t  he  presence  of  discrete  dispersal  delays. 
In  this  case,  large  values  of  k  give  rise  to  quasi -per  iodic  and  aperiodic  behavior,  but 
minimal  population  densities  remain  well  above  zero  for  both  prey  (Fig.  6-4 B)  and 
predator  (Fig.  6-4D). 


155 


Figure  (i-5:  Shape  of  the  Erlang  distribution  for  increasing  values  of  r.  I’lie  mean 
of  each  distribution  is  fixed  at  rav  —  2;  thus  b  —  c/2  for  each  curve.  Notice  that 
distribution  of  travel  times  is  narrower  for  larger  c. 

6.3.2  Distributed  travel  time 

When  the  movement  abilities  of  the  predators  differ,  or  t  he  vagaries  of  dispersal  affect 
individuals  differently,  individual  travel-times  form  some  distribution.  For  mathemat¬ 
ical  convenience  we  si  udy  a  ease  where  Ihe  delay  distribution  is  an  Erlang  distribution 

bc  sc~ 1  e-6s 

9(s)  =  SbAs )  =  jjj  •  (6.18) 

with  shape  parameter  c  and  scale  parameter  b  (Fig.  (i-5).  For  r  =  1  t he  distribution  is 
exponential,  and  the  dispersal  model  is  equivalent  to  one  that  includes  an  explicit  pool 
of  dispersers  with  constant  per  capita  emigration  and  immigration  rates,  a  la  Weisser 
&  Hassell  (1996)  and  Weisser  ct  al.  (1997).  For  c  >  l  the  mode  of  the  distribution, 
at  (c  —  1  )/b,  is  positive.  For  large  c,  the  mode  approaches  the  mean,  t„„  —  c/b,  and 
the  distribut  ion  resembles  a  delta,  function. 

For  this  special  family  of  distributions,  we  can  determine  the  local  stability  of  the 
equilibrium  point  (6.6)  by  linearizing  system  (6.5)  and  using  the  Laplace  transform 

f/M  =  9bAx)  =  )-b.  ,\c-  (6. 19) 

[x  +  by 
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Figure  6-6:  A-C)  Stability  diagrams  for  model  (6.5)  with  Erlang  distributed  dispersal- 
delays  computed  from  (6.20);  e  =  0.01,  /i  =  fip  =  1,k  =  1000.  A)  c  =  1.  B)  c  —  4,  C) 
c  =  64.  D)  Stability  diagram  for  the  Type  11  model  (6.17)  with  Erlang  distributed 
delay;  r  =  64.  Note  the  changing  scale  of  the  Tou-axis. 
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The  characteristic  equation  (6.14)  then  reduces  to  a  polynomial  of  degree  c  +  2: 


[(A  +  //  +  d)  (Xp*  —  sn*)  +  n*  (1  —  A)]  ^ —  d(Xp‘  —  en*) ,  (6.20} 

To  construct  the  stability  diagrams  in  Fig.  6-6.  we  found  the  roots  of  equation  (6.20) 
numerically  using  the  Matlab  function  rootsC). 

In  Figs.  6-6A-C  we  show  stability  diagrams  for  increasing  values  of  the  shape 
parameter  c,  with  the  other  parameters  fixed  at  levels  that  produce  an  unstable 
equilibrium  in  t  he  absence  of  dispersal  delays. 

The  area  in  the  parameter  space  where  the  equilibrium  point  is  stable  is  largest  for 
c  =  1  (Fig.  6-6A).  implying  that  the  stabilizing  effect  of  dispersal  delays  is  strongest, 
when  individual  travel  times  are  exponentially  distributed.  As  c  increases,  the  sta¬ 
bility  region  shrinks  and  its  borders  become  more  convoluted.  The  area  of  stability 
remains  large  even  in  the  limit  as  k  — >  oo  (Fig.  6-6D).  We  expect  that  as  c  becomes 
even  larger  the  stability  diagram  would  look  even  more  like’  the  discrete-delay  ease. 
Unfortunately,  for  r  much  larger  than  64.  the  characteristic  polynomial  (6.20)  is  ex¬ 
tremely  poorly  conditioned,  with  coefficients  differing  in  magnitude  by  hundreds  of 
orders;  we  have  been  unable  to  construct  a  stability  diagram  for  these  cases. 


6.4  Discussion 

We  have  shown  that  the  coexistence  equilibrium  point  of  the  single-patch  Rosen zweig- 
MacArthur  model  (6.2)  can  be  stabilized  when  predator  dispersal  includes  a  dispersal 
delay.  Stabilization  occurs  because  the  dispersal  delays  introduce  density  dependence 
into  the  dispersal  process  (Murdoch  et  al. ,  1992;  Neubert  ct  ai,  2002).  If  the  predator 
population  on  the  patch  is  abundant  compared  to  earlier  times,  emigration  from  the 
patch  will  exceed  immigration,  and  the  abundance  on  the  patch  decreases.  If,  on 
the  other  hand,  the  current  population  on  the  patch  is  small,  immigration  will  exceed 
emigrat ion,  thereby  increasing  the  population  size.  In  this  way,  population  oscillations 
are  reduced  and  species  abundances  eventually  reach  their  equilibrium  levels. 

The  stabilizing  effect  of  a  dispersal  delay  is  strongest  when  t  he  individual  travel 
times  are  exponentially  distributed,  as  they  are  in  models  that  include  a  pool  of 
dispersers.  The  stabilizing  effect  weakens  as  the  delay  distribution  becomes  more 
concentrated  around  its  mode.  In  the  weakest  case,  when  the  delay  distribution 
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is  a  delta  function  and  the  carrying  capacity  is  infinite,  the  stability  region  takes 
the  shape  of  an  archipelago  (Fig.  6-3B-D).  The  same  structure,  dubbed  “islands  of 
amplitude  death,"  has  been  observed  in  mathematical  studies  of  coupled  oscillators 
(Reddy  et  al,  1998,  1999).  These  studies  find  that  the  amplitude  of  two  coupled 
limit-cycle  oscillators  can  be  “quenched"  when  the  coupling  is  time  delayed. 

Dispersal  delays  are  less  effective  at  stabilizing  the  equilibrium  as  the  carrying 
capacity  of  the  prey  increases,  Nevertheless,  for  a  significant  set  of  parameter  values, 
the  model  with  dispersal  delays  has  a  stable  coexistence  equilibrium  even  for  an 
infinite  carrying  capacity  (see  Fig.  6-3  and  Fig.  6-6D).  Thus  dispersal  delays  alone 
are  capable  of  inducing  stability  in  the  face  of  a  destabilizing  Type  II  functional 
response. 

Even  when  the  equilibrium  is  unstable,  the  amplitude  of  the  predator-prey  oscil¬ 
lation  does  not  grow  with  increasing  carrying  capacity,  and  the  minimum  population 
densities  remain  well  above  zero  (Fig.  6-4B,  and  D).  In  tills  sense,  our  results  can  be 
added  to  those  of  Jansen  (1995)  (see  also  de  Roos  et  al.  1991;  Scheffer  V  de  Boer. 
1995;  Nisbet  et  al. ,  1998;  Jansen  k  de  Roos,  2000;  Jansen,  2001)  who  also  proposed 
dispersal  (without  delay)  as  a  potential  resolution  of  the  paradox  of  enrichment .  The 
st  abilizing  effect  of  dispersal  in  these  studies  is  weaker  than  it  is  in  our  model,  how¬ 
ever,  as  it.  only  produces  a  decrease  in  the  amplitude  of  the  limit  cycle,  rather  than 
stabilizing  the  equilibrium  point. 

Spatial  structure  is  by  no  means  the  only  factor  that  has  been  proposed  to  re¬ 
solve  the  paradox  (Abrams  k  Walters,  1996).  Other  factors  include  heterogeneity 
within  the  prey  population  and  complex  food  web  structure.  Enrichment  of  the  prey 
can  reduce  the  amplitude  of  population  cycles  when  prey  have  different  profitability 
(Genkai-Kato  k  Yamamura,  1999)  or  when  a  single  predator  attacks  two  prey  species, 
one  of  which  is  inedible  (Kretzschmar  et  al.  1993).  Enrichment  can  even  lead  to  sta¬ 
bility  in  systems  that  have  a  prey  refuge  (Abrams  k  Walters,  1996;  Gurney  k  Veitch, 
2000)  or  inducible  defences  in  prey  (Vos  et  al ,  2004).  Enhanced  system  persistence 
and  stability  in  intricate  food  webs  has  been  attributed  to  weak  trophic  interactions 
that  dampen  oscillations  between  consumers  and  resources  and  maintain  population 
densities  further  away  from  zero  (McCann  et  al. ,  1998). 

Our  analysis  has  several  limitations.  We  focussed  on  a  single  habitat  patch  from 
which  only  predators  dispersed.  Using  methods  outlined  Neubert  et  al  (2002),  one 
can  show  that  if  dispersal  delays  stabilize  the  single  patch  model  they  also  stabilize  a 
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spat  ially  homogeneous  equilibrium  of  a  model  with  an  arbitrary  number  of  identical 
patches.  Many  real  metapopulations,  however,  are  composed  of  numerous  patches 
that  differ  in  several  attributes,  in  particular,  the  distance  between  two  patches, 
and  therefore  the  distribution  of  dispersal  delays  between  them,  will  not  be  the  same 
for  all  pairs  of  patches.  Furthermore,  both  prey  and  predators  may  disperse.  Our 
analysis  docs  not  apply  to  these  more  complicated  scenarios. 

Finally,  we  note  that  our  results  may  depend  upon  the  exact  way  in  which  wc 
modeled  the  dispersal  process.  Another  approach  uses  so  called  “patch  occupancy 
models.”  which  keep  track  of  the  number  of  habitat  patches  that  are  in  various  states, 
e.g.,  empty,  or  occupied  by  prey,  or  occupied  by  predators.  In  contrast  to  our  results. 
Sabelis  ei  al.  (1991)  showed  that  while  the  addition  of  a  pool  of  dispersing  prey 
was  stabilizing  in  a  simple  patch  occupancy  model,  dispersing  predators  could  be 
destabilizing.  When  ii  comes  to  the  effects  of  dispersal  on  predator- prey  dynamics, 
the  details  of  how  dispersal  is  incorporated  appear  to  be  important  . 
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Chapter  7 


Future  directions:  Matrix 
population  models  for  epidemics 
and  demography 

7.1  Introduction 

Infectious  diseases  often  affect  host  individuals  of  different  (st)ages  in  a  different  way. 
For  example,  children  aged  6  59  months,  pregnant  women,  and  persons  aged  50  years 
or  more  are  at  much  higher  risk  for  influenza- related  complications  and  severe  disease 
than  persons  between  ages  5  and  50.  Different  stages  can  have  different  susceptibility 
to  a  disease  also  in  the  case  of  wildlife  infections.  In  case  of  the  phocine  distemper 
virus,  different  groups  have  different  behavior,  influencing  their  probability  of  getting 
infected. 

The  demographic  time-scale  is  usually  very  different  from  the  epidemic  time-scale, 
so  most  models  focus  either  on  demographic  or  epidemic  questions.  But,  in  order  to 
know  how  recurring  epidemics,  or  epidemics  that  have  long  infectious  period,  such 
as  HIV/AIDS,  affect  the  population,  the  model  should  incorporate  bot  h  realistic  de¬ 
mographic  and  epidemic  detail,  One  way  to  incorporate  age  structure  into  epidemic 
models  has  been  by  using  integro-differential  equations  (e.g.,  Diekmann  &  Heester- 
beek,  2090;  Hethcote,  2000;  Dietz  k  Heesterbeek,  2002;  Thieme,  2003).  Another  way 
is  to  use  matrix  models.  Apart  the  applications  in  the  analysis  of  age-prevalence  data 
(Saporu,  1990.  1996),  and  description  of  the  contact  and  transmission  processes  (e.  g., 
Pugliese,  1991;  Keeling  &  Grenfell,  1997).  the  matrix  approach  hasn’t  commonly  been 
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used  in  epidemiology. 


I  want  to  construct  a  matrix  model  that  accounts  for  both  epidemic  and  demo¬ 
graphic  detail,  and  look  at  its  dynamics. 


7.2  Epidemic  model  with  demography 


To  allow  for  demographic  detail  in  each  of  the  epidemic  compartment  1  will  follow  the 
approach  formulated  by  Hunter  Caswell  (2005)  to  model  spatial  matrix  population 
models.  The  authors  constructed  matrix  models  from  a  manageable  block-diagonal 
formulation  of  the  dispersal  and  demographic  processes,  using  a  special  permutation 
matrix  called  a  vec-permutatiou  matrix. 

Even  though  the  spatial  spread  of  infectious  diseases  is  an  important  area  of  re¬ 
search  today,  here  I  am  not  interested  in  spatial  aspect  of  an  epidemic,  Instead,  1  want 
to  apply  the  vec- permutation  approach  to  st  udy  epidemic  processes  in  a  demographic 
sett  ing. 

Let  tin1  host  population  be  divided  into  s  stages  and  c  epidemic  categories.  The 
total  number  of  population  compartments  is  then  s  x  e.  The  state  of  the  epidemic  in 
tlu'  population  at  time  /  can  then  be  described  by  the  matrix 
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where  rijj(f)  is  the  number  of  individuals  in  stage  i  and  in  epidemic  category  j  at  time 
1.  Row  i  (u,.).  for  example,  has  individuals  of  the  same  stage  but  in  different  epidemic 
categories,  where  column  j  (?/.,  )  gives  individuals  in  the  same  epidemic  category,  but. 
in  different  stages.  Forming  a  population  vector  by  stacking  the  rows  of  this  matrix, 
will  give  us  all  the  individuals  in  stage  1,  followed  by  all  the  individuals  in  stage  2, 


et  c. 
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Accordingly,  we  can  stack  the  population  vector  according  to  the  epidemic  categories, 
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The  population  vector  n  can  be  easily  organized  into  epidemic  or  demographic 
stages  by  the  use  of  the  vec  operator,  vec(>),  that  stacks  the  columns  of  a  matrix  on 
top  of  each  other, 


"stages  =  vec(iVT), 


nepidemic  vec  (N). 


The  vectors  (7.4)  and  (7.5)  are  related  by  a  special  matrix  called  the  vec-permutatian 
matrix  P  so  that 


vec(A,T)  =  P  vec(N). 

Since  P  is  a  permutation  matrix,  it  holds  t  hat 


PT  =  P"1. 


(7.7) 


Imagine  a  case  where  a  population  reproduces  at  the  beginning  of  the  projection 
interval.  As  we  want  to  keep  track  of  both  the  demographic  stages  and  epidemic 
categories  of  the  individuals,  upon  birth  newly  born  individuals  are  assigned  to  their 
epidemic  compartments.  An  infection  occurs  at  the  end  of  the  projection  interval, 
and  let  the  matrix  A[n(i)]  describe  the  epidemic  transitions  that  occur  during  this 
outbreak.  Let  the  matrix  K  describe  reproductive  events,  and  the  matrix  M  move  the 
newborn  individuals  to  appropriate  epidemic  categories.  Using  the  ver-pernmtation 
configuration  and  a  population  vector  organized  into  epidemic  categories,  we  can 
describe  the  transitions  that  occur  in  one  projection  interval  by 
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Alternatively,  an  outbreak  of  infection  may  occur  at  the  beginning  of  the  projec¬ 
tion  interval,  followed  by  the  reproduction,  and  ‘movement.’.  In  this  cast;,  it  is  more 


convenient  to  organize  the  population  according  to  its  demographic  compartments, 


(<+!)  =  MPRPTA[n(f)] 


(0- 


(7.9) 


7.2.1  A  (not  so)  simple  example 

To  see  how  to  formulate  epidemic  and  demographic  transition  matrices,  consider  a 
population  that  consists  of  only  juveniles  and  adults  2  demographic  stages,  s  =  2. 
Further  assume  that  the  disease  that  invades  this  population  follows  the  standard 
SIR -type  dynamics,  so  at.  any  given  time  an  individual  can  be  either  susceptible  (5), 
infectious  (/),  or  recovered  (/?)  3  epidemic  categories,  c  =  3.  In  this  example,  the 

total  number  of  compartments  in  the  population  is  sxc=  0.  The  population  matrix 
can  be  described  as 


SIR 


(7.10) 


Epidemic  part 

Let  the  t  ransmission  rate  j3  vary  between  groups,  so  that  St ,  is  the  transmission  rate 
between  a  susceptible  individual  in  demographic  stage  i  and  an  infectious  individual 
in  demographic  stage  j.  The  probability  that  a  suscept  ible  in  stage  1  has  no  infectious 
contacts  with  infectives  in  stage  1  during  the  time  interval  (/,/  +  !)  is  exp  (— 
Accordingly,  the  probability  that  a  susceptible  in  stage  1  has  no  contacts  with  infec¬ 
tives  in  stage  2  is  exp  (—  0i2ri22(t)).  The  total  force  of  infection  A|[n(/)]  for  stage  1 
(juveniles)  at  time  t  is  then 


(7.11) 


After  getting  infected,  a  individual  in  stage  i  remains  infectious  for  an  average  dura¬ 
tion  of  1/7;,  after  which  it  recovers  with  the  probability  r,.  We  can  summarize  the 
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epidemic'  transition  with  the  matrix  A(n(f)] 


A[n(/)]  = 


^1-Ai[n(()]  0  0 

\ 

Ai[n(0]  l-7i  0 

0 

0  ri7i  1 

1-Ai|n(0]  0  0 

0 

A2[n(t)]  I-72  0 

0  r2 72  1  j 

and  we  can  project  the  epidemic'  via 

n(f  +  1)  =  A[n{f)]n(0- 


(7.12) 


(7.13) 


When  we  don’t  account  for  the  demographic  detail,  an  outbreak  of  disease  quickly 
grows  into  an  epidemic,  and  finally  disappears  from  the  population  (see  Figure  7-4). 
Without  the  reproduction,  there  is  no  influx  of  susceptibles,  so  hosts  are  quickly 
exhausted  and  the  disease  disappears  from  the  population  without  the  possibility  of 
reaching  an  endemic  equilibrium. 


Demographic  part 

The  population  at  any  time-step  consists  of  juveniles  (nj.)  and  adults  {n^. )  that,  can 
either  be  susceptible  (n.i),  infectious  (n.a).  or  recovered  (n.;f).  Individuals  in  the 
stage  i  and  epidemic  category  j  suffer  mortality  m from  natural,  non-disease  related 
causes.  Juveniles  in  the  epidemic  category  i  survive  and  grow  to  adults  (in  the  same 
epidemic  category)  with  the  probability  (/*.  Adults  in  the  epidemic  category  i  have  the 
per-capita  fertility  /,.  and  they  produce  newborns  that  arc  temporarily  in  three  new 
demographic  stages,  n:(.  (see  Figure  7-2  for  an  illustration).  Tin'  temporary  stage 
7>;ti  consists  of  newborns  produced  by  susceptibles,  n 32  are  newborns  produced  by 
infect  eds.  anti  recovered  individuals  give  birth  to  n:«. 

If  we  summarize  the  transitions  of  epidemic  category  i  with  the  matrix  R, . 


R,  - 


(  (l-ft)(l  ~  rnu)  0  ^ 
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The  reproduction,  survival,  and  growth  can  lie  written  using  the  block-diagonal  form 
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Figure  7-1:  Dynamic  of  the  epidemic  in  a  structured  population  without  demography. 
Parameter  values:  =  0.005.  71  =  7 2  =  1/14,  ri  =  r2  =  0.5.  initial  conditions: 

nstages  =  (  100  0  0  50  1  0  ) 


susceptibles  infcctives  recovered 


Figure  7-2:  Reproduction  matrix  R  accounts  for  the  reproduction,  survival  and 
growth.  New  individuals  temporarily  show  up  in  extra  states  (n:t.)  before  they  are 
assigned  to  their  epidemic  categories. 
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juveniles 


newborns 


adults 


vertical  inherited 

transmission  (v)  immunity  (h) 


Figure  7-3:  After  reproduction,  newborn  individuals  are  assigned  to  their  epidemic 
categories  by  the  matrix  ML  If  there  is  no  vertical  transmission  or  inherited  immunity 
(v  =  0,  //  =  0)  then  all  newborn  individuals  will  be  susceptible. 


;LS 


( tin  ^ 
nat 

nii 

nn 

nn 

«13 

n.23 

«33  / 


n„  \ 


Ri 

0 

0  \ 

0 

Ri 

0 

0 

0 

R;<  / 

j 

»23  / 


(7.15) 


where  K,  is  a  9  x  (5  matrix  that  consists  of  3  x  2  blocks  R,  on  the  diagonal  and  zeros 
elsewhere. 


Newly  born  individuals  do  not  have  to  be  in  the  same  epidemic  category  category 
as  their  parents.  If  there  is  no  vertical  transmission  of  a  disease,  new  individuals  born 
to  infectious  parents  will  be  susceptible.  In  case  that  newborns  don’t  have  maternal 
antibodies  (i.  e.,  inherited  immunity)  against  an  infection,  even  though  they  are  born 
to  immune  (recovered)  parents,  they  will  be  susceptible.  Figure  7-3  summarizes  those 
possibilities. 

Assignment  of  newborn  individuals  to  epidemic  categories  can  be  summarized  by 
the  ‘movement’  matrix  M  and  the  following  equation. 
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population  projection  interval 

Figure  7-4:  Without  the  epidemic,  the  structured  population  grows  exponentially. 
Parameter  values:  g\  =  0.3.  m\.  —  0.2.  m2-  —  0.1.  /,  =  0.3.  initial  conditions: 

■■epidemic  =  (  1«)  50  0  0  0  0  )'. 
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Iii  the  absence  of  infections,  the  dynamics  of  the  population  is  given  by 

n(/.  +  l)=P6TMP()Mn(/.)  (7.17) 

where  and  l\ t  are  permutation  matrices  described  in  (7.0),  Without  the  epidemic, 
the  population  grows  exponentially,  without  any  infectious  individuals.  The  pop¬ 
ulation  growth  rate  is  given  by  the  dominant  eigenvalue  A|  of  the  square  matrix 
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Figure  7-5:  Endemic  equilibrium  of  model  (7.18)  that  includes  both  demographic 
and  epidemic  detail.  Parameter  values:  =  0.3.  m|.  =  0.2,  m2.  =  0.1.  /,  =  0.3. 

Jij  =  0,005.  7)  =  72  =  1/14,  77  =  r 2  =  0.5,  h  =  0,  v  =  0.  initial  conditions: 

^epidemic  =  (  50  50  1  0  0  0  )T. 

B  =  Pf[  M  Pt)  R. 

Model  with  demography  and  epidemics 

Combining  demographic  and  epidemic  detail  can  be  done  in  several  ways.  For  exam¬ 
ple.  the  reproduction  can  be  followed  by  an  outbreak  of  a  disease  as  in  (7.8).  or  there 
can  be  an  outbreak  at  the  beginning  of  the  population  projection  interval,  followed 
by  the  reproduction  as  in  (7.9). 

Consider  a  case  where  an  epidemic  occurs  at  the  end  of  the  population  projection 
interval. 

n{t  +  l)  =  PfI  A[n(i+Mi)]-A[n((  +  A/)jA[n(t)]  MP„Rn(/)  (7.18) 

The  epidemic  introduces  nonlinearity  into  the  model  that  brings  this,  otherwise 
exponentially  growing,  population  to  an  equilibrium.  The  demographic  part,  on  the 
other  hand,  introduces  an  influx  of  susceptible®  into  the  population,  which  allows 
the  disease  to  become  endemic.  Combining  a  standard  SIR  dynamic  with  simple 
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demographic  model  allows  a  new  type  of  behavior  -  endemic  equilibrium  -  which  was 
not  possible  in  models  (7  13)  and  (7.17). 

7.3  Future  analyses 

Combining  epidemic  and  demographic  detail  into  a  single  model,  gives  rise  to  the 
dynamics  not  present  in  the  building  blocks  of  this  model.  In  the  future  1  want  to 
explore  the  dynamics  of  this  model  in  more  detail  and  look  at  the  endemic  equilibria 
and  the  stability  of  this  system.  What  epidemic  and  what  demographic  parameters 
can  drive  this  system  to  instability?  I'm  interested  to  see  under  what  conditions  can 
a  disease  persist  in  this  model,  and  what  factors  lead  to  disease  extinction. 
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Appendix  A 

Proof  that  system  (5.19)  has  no 
real  solution 


To  rule  out  imaginary  roots  for  the  characteristic  equation  (5. 14)  wo  prove  that  system 
(5.19)  has  no  real  solution.  To  complete  the  proof  we  use  a  theorem  by  Hardy  et  al. 
(1952).  To  use  the  this  theorem  we  first  need  the  following  definitions. 


Definition  1  Max f.  the  ’ effective  upper  bound'  of  /,  is  defined  to  be  the  largest  £ 
which  has  the  following  properly:  if  <  >  0,  there  is  a  set  r(r )  of  positive  measure  in 
which  f  >  £  —  f  .  If  there  is  no  such  £,  we  write  Max /  —  oc.  For  funct  ions  continuous 
on  a  closed  interval.  Max  f  is  the  ordinary  maximum. 


Definition  2  The  mean  of  f  with  respect  to  the  weight  function  0  on  a  nuns  arable 
set  E  is  defined  as 


W(/)  - 


0(*)  /(«)  ‘l* 

jEcp{s)ds 


(Al) 


With  these  definitions,  we  can  now  state  the  theorem. 


Theorem  1  (Hardy  et  al.  1952,  Theorem  183)  Let  the  measurable  function  f 
be  finite  almost  everywhere  on  a  measurable  set  E  and  non-negative.  Let  the  measur¬ 
able  function  0  be  finite  and  positive  everywhere  in  E,  and  integrable  over  E.  Then, 
if  U(  f)  is  finite  and  positive, 

U{f)  <  Max/,  (A2) 

unless  f  =  C  ( C  a  constant)  almost  everywhere. 


We  now  prove  the  following  lemma. 


Lemma  1  If  dp.  fi.  n*  and  p*  arc.  positive,  pv  is  nonnegative,  and  the  probability 
density  function  gp(s)  is  finite  and  has  support  on  a  measurable  set,  then  system 
(5.1(J)  has  no  real  solution  y. 


Proof.  One  solution  of  equation  (5,19b)  is  y  =  0.  But  y  =  0  is  not  a  solution  of 
(5.19a),  so  y  ^  0.  Dividing  equation  (5.19b)  by  dpygp(p,p)  gives 


which  implies 


£  <lA*)  e  w  eos(ys)  (Is 
£  9P(s)  e-*”*  ds 

.£jh£)££  |cos(,^)|r/,s 

£  9p(s)  e-11”*  ds 


(All) 


(A4) 


But  Theorem  1.  with  d{s)  =  gp(s)e~>1,’a .  f(s)  —  |eos(y.s)|,  E  —  {s  >  0  :  gp{&)  >  9} 
and  y  ^  0.  implies  that  the  left  hand  side  of  equation  ( A4 }  is  less  than  one,  Thus 
system  (5.19)  does  not  have  a  solution.  □ 
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